General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 





CML 82-5 


NASA C r : 68057 


WAVE PROPAGATION IN GRAPHITE/EPOXY 
LAMINATES DUE TO IMPACT 

by 

T.M. Tan and C.T. Sun 


December, 1982 


1 Htpnrt No 2. Oanornmoni A cc ession No 3. Recipients Catalog No 

NAS* CR- 16805 7 


4 Tula and SuMnii 

S. Papon Dim 

WAVE PROPAGATION IN GRAPHITE/EPOXY LAMINATES DUE TO IMPACT 

December 1982 


6 Performing Organization Copt 

7 Author In 

•. Performing Organization Raport No 

T. M. Tan and C. T. Sun 



10. Work Ural No 

i Performing Organization Name r Address 


Purdue University 


School of Aeronautics and Astronautics 

IV Contract or Grant No 

Nest Lafayette, IN 47907 

NSG 3185 


13. Type or Raport and Par tod Cowar ad 

13. Sponsoring A fancy Nam* and Address 

Interim Report 

National Aeronautics & Space Administration 


Washington, OC 20546 

14. Sponsoring A fancy Coda 

1 is Supplementary Notes Project Monitor: C. C. Chamis, Structures (Mechanical Technologies Dlv. 

NASA Lewis Research Center, M.S. 

49-6 

21000 Brookpark Road 


Cleveland, OH 44135 


16 Abstract 


The low velocity impact response of graphite/epoxy laminates Is investigated theoretically and 
experimental ly. A 9-node Isoparametric finite element In conjunction with an empirical contact 
law was used for the theoretical Investigation. Flat laminates subjected to pendulum Impact were 
used for the experimental Investigation. Theoretical results are In good agreement with strain 
gage experimental data. The collective results of the Investigation indicate that the theoretical 
procedure describes the impact response of the laminate up to about 150 In/sec. impact velocity. 


ORIGINAL PAGE 1$ 
OF POOR QUALITY 


17 Kay Words I Suggested by Authorial I 

Impact fiber composites, laminates contact law, 
finite element, stress waves 

11. Distribution Statement 

Unclassified, Unlimited 

19 Security Oattil (of this raportl 

Unclassi fled 



20. Security Classit (of this papal 
Unlimited 

21 No ol Pages 

130 

22. Pt'ca* 


’ For sale by the Na' onal Technical Information Service Spnngneifl Virginia 22161 


NASA-C-168 (Rpv 10-75) 
























m 


ORIGINAL PAGE IS 
OF POOR QUALITY 


TABLE OF CONTENTS 


Pag© 


TABLE OF CONTENTS Ill 

LIST OF TABLES v 


LIST OF FIGURES Vi 

LIST OF SYMBOLS lx 

CHAPTER i - INTRODUCTION . 1 

CHAPTER 2 - STRESS WAVE IN A LAMINATED PLATE 5 

2.1 Laminate Theory with Transverse shear effects 6 

2.1.1 Lamina Constitutive Equations 6 

2.1.2 Plate Strain-Displacement Relations . 9 

2, .1.3 Stress-Resultants and Laminate Constitutive 

Equations 12 

2.1.4 Plate Equations of Motion 16 

2 . 2 Propagat l on of Harmon I c Waves . 18 

2.3 Propagation of Wave Front 23 

2.3.1 Kinematic Conditions of Compatibility on the 

Wave Front 25 

2.3.2 Dynamical Conditions on the Wave Front ...... 27 

2.3.3 Propagation Velocity of the Wave Front 32 

2.3.4 Wave Surface and Ray 35 

CHAPTER 3 - STATICAL INDENTATION LAWS 48 

3.1 Specimens and Experimental Procedure 52 

3.2 Experimental Results 53 

3.2.1 Loading Curves 53 

3.2.2 Unloading Curves 57 

3.2.3 Reloading Curves 71 


preceding page blank not filmed 



! v 

3.3 Discussion . . 71 

CHAPTER 4 - IMPACT EXPERIMENTS 80 

4. 1 Experimental Procedure 81 

4.2 Calibration of Impact-Force Transducer 84 

4.3 Finite El ement Ana lysis 94 

4.3.1 Plate Finite El ement 04 

4.3.2 Modeling of Projectile 97 

44 Results and Discussion 100 

CHAPTER 5 - SUMMARY AND CONCLUSION 113 

LIST OF REFERENCES 116 


APPENDIX: COMPUTER PROGRAM AND USER INSTRUCTIONS 


119 


V 


LIST OF TABLES 

Table 

Page 


3.1 Contact coefficient l< of loading law F * 60 

4.1 Specifications for Model 200A05 Impact-Force 

Transducer 85 



vl 


LIST OF FIGURES 


Figure 


Page 


2.1 Lamina reference axes and laminate reference 

axes 8 

2.2 Laminate displacement components for a 

cross-section perpend Icul ar to the y-axls. , 10 

2.3 Stress-resul tants and geometry of a typical 

N- 1 ayer 1 am I nate . . 14 

2.4 Dispersion curves for plane harmonic waves 
propagating In the 0°~ 45°- and 90°~ 

directions. 22 

2.5 Frequency curves for flexural waves 
propagating In the 0°~ 45°- and 90 u - 

d I rect I ons 24 

2.6 A deformed volume V divided by a travelling 

surface Q. , 29 

2.7 Normal velocities of In-plane wave fronts... 36 

2.8 Normal velocities o^ flexural wave fronts..,, 37 

2.9 Wave front positions at different times and 

rays for In-plane extenslonal mode 43 

2.10 Wave front positions at different times and 

rays for In-plane shear mode 44 

2.11 Wave front positions at different times and 

rays for bending mode 45 

2.12 Wave front positions at different times and 

rays for twisting mode... 46 

3.1 Schematlcal diagram for the Indentation test 

set-up 54 

it 

3.2 Loading curve of [O 0 /45 0 /0°/-45 0 /O 0 ] 2S 

specimens with 0.5 Inch Indenter (n E =3/2) 55 



vl I 


3.3 Loading curve of [90 0 /45°/90 0 /~45 0 /90 0 ] ?s . 
specimens with 0,5 Inch Indenter (n-3/2). 


* * * » * » • * * 


56 


3.4 Loading curve of [0°/45 0 /0°/-45 0 /0 0 ] as 

specimens with 0.75 Inch Indenter (n*3/2) 58 


3.5 Loading curve of [90°/45 0 /90°/-45 0 /90 0 3 as 

specimens with 0.75 Inch Indenter (n*3/2) , 59 


3.6 Relation between permanent Indentation and 

maximum Indentation . 61 


3.7 Unloading curves of [O°/45 o /O o /-45 0 /O°] as 

specimens with 0.5 Inch Indenter (q»2.2) , 63 

3.8 Unloading curves of [90 c /45 0 /90°/“45 0 /90°3 2S 

specimens with 0.5 Inch Indenter (q-2,2) 64 

3.9 Unloading curves of [0°/45 o /0 o /-45 o /0°] 2S 

specimens with 0.75 Inch Indenter (q«1.8) 65 


3.10 Unloading curves of [90°/45 0 /90 0 /-45°/90 0 ] 2S 

specimens with 0,75 Inch Indenter ( q*1 =8) = ,.,,..-,; 66 


3.11 Unloading curves of [0 0 /45 0 /G 0 /~45°/0 0 ] 2S 

specimens with 0.5 Inch Indenter (q=»2.5) .......... 67 

3.12 Unloading curves of [90°/45 0 /90 0 /“45°/90 0 ] 2S 

specimens with 0.5 Inch Indenter (q**2,5) . , . . 68 

3.13 Unloading curves of [0°/45 0 /0 0 /~45 0 /0 0 ] 2S . 

specimens with 0.75 Inch Indenter (q-2.0) 69 


3.14 Unloading curves of [90°/45 o /90 o /-45°/90 o ] 2S 

specimens with 0.75 Inch Indenter (q=2.0) 70 

3.15 Reloading curve of [O 0 /45 0 /0°/~45 0 /0 0 ] 2S 

specimen with 0.5 Inch Indenter (p-1.5) 72 


3.16 Reloading curve of [90 0 /45°/90 0 /-45 0 /90 0 ] 2S 

specimen with 0,5 Inch Indenter (p“1.5)... 73 

3.17 Reloading curve of [O 0 /45 0 /O 0 /~45 0 /O 0 ] 2S 

specimen with 0.75 Inch Indenter (p=1 .5) 74 

3.18 Reloading curve of [9O 0 /45 o /9O°/-45 o /9O o 3 2S 

specimen with 0.75 Inch Indenter (p=1 .5) 75 

3.19 Unloading rigidity s as function of maximum 

l ndentat I on. 78 


4.1 Laminate dimension and strain gage locations 


82 


vtll 


4.2 Graphical Illustration of Impact projectile,,,,.,, 82 

4.3 Schematlcal diagram for the Impact 

experimental set-up. ,**»»*»*»•»»♦,,»»»•«»*,»*«•»*» 83 

4.4 Experimental set-up for the calibration of 

Impact-force transducer . . . 86 

4.5 Typical output voltages from transducer and 

accel erometer 88 

4.6 Relation between V F and V a , . 88 

4.7 Assumed exponential Impulsive loading and 
the response history at the midpoint of the 

rod 90 

4.8 Accelerations of rod for assumed exponential 

Impulsive loading 92 

4.9 Assumed s I ne-funct 1 on Impulsive loading and 
the response history at the midpoint of the 

rod 93 

4.10 9-node 1 soparametr lc plate element., 95 

4.11 Finite element mesh for laminated plate and 

project! le 101 

4.12 Strain response history at gage No.1 103 

4.13 Strain response history at gage No. 2 104 

4.14 Strain response history at gage No. 3 105 

4.15 Strain response history at gage No. 4 106 

4.16 Strain response history at gage No. 5....,., 107 

4.17 Strain response history at gage No. 6 108 

4.18 Transducer response and contact force 
histories from experimental and finite 

element results 109 

4.19 Transducer response histories from 
experimental and finite element results up 

to 800 microseconds Ill 

4.20 Deformed configurations of laminated plate 

after impact 112 



lx 


n 


|H- 


LIST OF SYMBOLS 


A 

A i j . B t j , D | j 
Eg 

Ei 

E-2 

P 

F m 

G 

[K p 3 , EK r ] 

[M p j EM r ] 

M 

N 

{Pp}, { P r > 

Q 

Q.j 

Qi j 

R s 

s, 

V F 

V a 


Cross-sectional area of the projectile 

Laminate stiffnesses 

Young's modulus of the steel Indenter 

Young's modulus of laminar In the fiber 
direction 

Young's modulus of laminar In the transverse 
direction 

Contact force 

Maximum contact force 

Shear modulus 

Stiffness matrices 

Mass matrices 

Stress couples of laminate 
Stress resultants of laminate 
Assembled global load vectors 
Transverse shear force of laminate 
Reduced stiffnesses 
Transformed reduced stiffnesses 
Radius of steel Indenter 
Shape functions of plate element 
Output voltage of the force transducer 
Output voltage of the accelerometer 


V ?* 

■ ’ i 




L L 


H h 

it i 


v 


O XV 





U it 




■ 4 } ; 


»r v 

» I! 

- IS 


I f 


U 



X 


a 

c 

C« 

Cp 

Cn 

f| 

Cf 3 

h 

k 

k 

ki 

Ek p 3 « Ek r 3 
[m P .l, [m r 3 
n 

n i 

P 

Pi 

{ Pp I a > EPr^o 

q 

C q p ) # { q r ) 

E q P 1 Q t f qr 3 o 

s 

t 

t* 

U, V, w 
u°, v°, w° 


Acceleration 
Phase velocity 

Sensitivity of the accelerometer 

Sensitivity of the Impact-force transducer 

Normal velocity of wave front 

Shape functions of rod element 

Discontinuity of f across *Bve front surface 

Laminate thickness 

Wave number 

Contact coefficient 

Reloading rigidity 

Element stiffness matrices 

Element mass matrices 

Power l ndex of I oad I ng law 

Unit normal on the wave front 

Power Index of reloading law 

Slowness vector 

Element load vectors 

Power Index of unloading law 

Assembled global displacement vectors 

Element displacement vectors 

Unloading rigidity 

Time 

Non-d I mens lonal time 
Displacement components of laminate 
MI dp lane displacement components 



xl 


x, y, z Laminate coordinate system 

x, , x 2 , K$ Laminar coordinate system 


n 

a 

«o 

01m 

r > 

V 

e 

^x » i ^xy 

K 

v 


Vs 

$. 7? 

P 

O’ 

r 

<Px » ^y 
00 


Wave front surface 
Indentation depth 
Permanent Indentation 
Maximum Indentation 
Critical Indentations 
Shearing strain 
Normal strain 
Rotation gradients 
Wave length 
Poisson's ratio 

Poisson's ratio of the steel Indenter 
Normalized local coordinates of plate element 
Mass density of laminate 
Normal stress 
Shearing stress 

Rotations of cross-sections of laminate 
Frequency 


1 

ORIGINAL PAGli VJ 
OF POOR QUALITY 


CHAPTER 1 
INTRODUCTION 

Advanced f Iber-relnforced composite materials such as 
boron/epoxy and graph I te/epoxy have been successfully 
employed as structural materials In aircrafts,, miss I les and 
space vehicles In recent years, and the performance of these 
composites has shown their superiority over metals In 
applications requiring high strength, high stiffness as well 
as low we I ght . The advantages of these compos l tes , however , 
are overshadowed by their relatively poor resistance to the 
Impact loadings, which has prevented the application of 
these materials to turbine fan bladings. Many other reports 
dealing with the responses of advanced composites to various 
types of Impact have further Increased the need for a better 
understanding of the problem so that the survivability of 
these composites can be Improved. 

It Is obvious that Impact produces damage and 
consequently reduces the strength of composite materials. 
The damage modes usuilly Include local permanent 
deformations, breakage of fibers, del ami nat Ions , etc.. 
While the cause of these damages are still unknown and may 
not be simple In nature, In general, the Impact of a soft 
object could give a longer contact duration, and the dynamic 
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response pf the whole structure Is of Importance, The hard 
object Impact usually gives a short contact time and results 
In the Initial transmlsson of Impact energy Into a local 
region of the structure. This Initial energy will propagate 
Into the rest of the structure In the form of stress waves, 
Far field damage away from the Impact area could result from 
the reflection of stress waves. It Is generally agreed that 
the cause of the sudden failure must be examined from the 
point of transient wave propagation phenomena, 

Flexural waves Induced by dynamic loads In laminated 
composites are more complicated than those In homogeneous 
and isotropic plates due to the anisotropic and 
nonhomogeneous properties In the laminate. Moreover, 
because of the low transverse shear modulus In fiber 
composites, the effect of transverse shear deformation 
becomes significant and should be considered In the 
formulat Ion. In Chapter 2, the laminate theory which 
Includes the transverse shear deformation effect Is 
reviewed, and harmonic waves In a graphite/epoxy laminated 
plate are studied. The propagation of wave front which, for 
a given time after Impact, bound the stressed region 
surrounding the Impact point, Is also Investigated. 

A survey of wave propagation and Impact In composite 
materials has been given by Moon [1]. Many analytical [2~ 
5], numerical [6—73 and experimental [ 8—1 01 methods have 
been employed to study the transient Impact problems. The 
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respone of a laminated plate can be analyzed using these 
methods provided the applied load history Is prescribed. 
However If the dynamic load results from an Impact of an 
object on the laminated plate, then the resulting contact 
force must be determined first. An accurate account of the 
contact behavior becomes the most Important step In 
analyzing the Impact response problems. 


A classical contact law between two elastic spheres was 
derived by Hertz [11]. When letting the radius of one of 
the spheres go to Infinity, one obtains the contact law 
between an elastic sphere and an elastic half-space. Many 
authors have used the Hertzian contact law for the study of 
Impact on metals and composites. [12-13], Recently, Yang and 
Sun [14] performed statical Indentation tests on graphite/ 
epoxy composite laminates using spherical steel indenters of 
different sizes and found that the Hertzian law of contact 
was not adequate. In particular, they found that 
significant permanent Indentations existed and that the 
unloading paths were very different from the loading path. 
Noting that energy dissipation takes place during the 
process of impact, Yang and Sun [14] suggested that this 
energy is responsible for the local damage of the target 
materials. The unloading curves and permanent Indentations 
obtained from the statical indentation tests may provide a 
useful information in estimating the amount of damage due to 
impact since this energy is simply the area enclosed by the 
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loading-unloading curves. In this study, similar statical 
Indentation tests were conducted and the results are 
presented In Chapter 3. 

Wang El 5] has performed a number of Impact tests on 
graph 1 te/epoxy laminated beams and plates. It was shown 
that the strain responses calculated using' finite element 
method and the statically determined contact laws from [143 
agreed with the experimental measurements quite well. This 
Indicates that the statical Indentation law Is reasonably 
adequate In the dynamical Impact analysis. It was also 
suggested that the contact force should be measured 
experimentally to provide an additional basis for comparison 
with the finite element solution which could allow further 
evaluation the applicability of the contact laws In Impact 
analysis. Chapter 4 describes an Impact experiment on 
graphite/epoxy laminated plate using an Impact-force 
transducer with a spherical steel cap as the Impactor. The 
contact force history and strain responses at various points 
on the plate were measured by means of the transducer and 
surface strain gages, respectively, and were compared with 
the predictions of finite element analysis using the 
statically determined contact law. 

Chapter 5 summarizes the results obtained In Chapter 2, 3 
and 4. 
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CHAPTER 2 

STRESS WAVE IN A LAMINATED PLATE 

A laminated plat© theory which Includes the effects of 
transverse shear deformation and rotatory Inertia was 
developed by Yang, Norris and Stavsky [16] in a way 
suggested by Mtndlin [17] for homogeneous Isotropic plates, 
It was shown that the frequency curves for the propagation 
of harmonic waves In an Infinite two- layer Isotropic plate 
In plane strain agreed with the predictions of the exact 
solution obtained from theory of elasticity very well, A 
similar laminated plate theory was developed by Whitney and 
Pagano [18] and was employed tn the study of static bending 
and vibration for anti symmetr lc angle^ply composite plates 
with particular layer properties. It was found that the 
effect of shear deformation can be quite significant for 
composite plates with span-to-depth ratio as high as 20. 
Good agreement was also observed In numerical results for 
plate bonding as comparing with exact solutions of 
elasticity. In this study, the laminate theory developed by 
Whitney and Pagano was used for Its simplicity yet quite 
satisfactory in describing the harmonic wave propagation 
[193, 
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2.1 Laminate Theory with Transverse Shear Effects 

2.1.1 Lamina Constitutive Equations 

A laminated plate of constant thickness h consists of a 
number of thin lamlnas of unidirectional ly f tber-relnforced 
composite perfectly bonded together. Each lamina, whose 
fiber may orient In any arbitrary direction, can be regarded 
as a homogeneous orthotropic solid. Consider a typical k~th 
lamina. A coordinate system (x 1 , x a , x a ) Is chosen In such 
a way that the Xf-x a plane coincides with the midplane of 
lamina, and x, and x a axes are parallel and perpend 1 cul ar to 
the fiber direction, respectively. If a state of plane 
stress parallel to the x t “X a plane Is assumed, then the In- 
plane stress-strain relations are given by 


k k 

r a. a a % 


*11 


*Qi i Qia 0 ' 


6 i i 

(T 2 2 

. fag 

Qi a Qa a 0 

• 

e a a 

y 1 2 . 


0 0 Qee. 


Tt a. 


(2-1 ) 


The transverse shear stress-strain relations are given by 



(2-2) 


In which 
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v 21 E 1 /( 1 *-J/ 12 v 21 ) 


(2-3) 


are the so-called reduced stiffnesses, where E, G and v are 
Young’s modulus, shear modulus and Poisson's ratio, 
respectively, and subscripts 1 and 2 denote the directions 
parallel to x., and x 2 axes, respectively. 

The coordinate system for an arbitrarily oriented lamina 
does not, In general, coincide with the reference axes 
(x,y,z) of laminated plate (see Figure 2.1). Using the 
coordinate transformat I on laws for stress and strain, we 
obtain the stress-strain relations In laminate reference 
system as 


k k 



In which <2,j are given by 


(2-4) 


§11 = Qi im 4 +2(Q, 2 +2Q SG )m 2 n 2 +Q 22 n 4 



8 


0 RlG\W&J- 

Of pOOvt 


i iVV 


QUPvtn 



I 


z,x 3 


X 


(X, yX^Xj) — Lamina Reference Axes 
(X y Y* Z ) — Laminate Reference Axes 


Figure 2.1 Lamina reference axes and laminate reference 
axes 
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$22 ■ Qi in 4 +2(Q 1 a +2Q 66 )m 2 n 2 +Q 22 m 4 
** (Q, i+Q 2 2~4Q 6G )m 2 n a +Q 1 a (m 4 +n 4 ) 

■ (Qi i ~Qi a-SOge )m 3 n+(Q 1 2“Q22 + 2Q GG )mn a 
Q 2 e * (Qi i“Qi 2 ~ 2 Q G6 )mn 3 +(Q 1 2“Qa2 + 2Q 66 )m 3 n (2—5) 

^66 * (Qi 1+Q22-2Q, 2 - 2 Q GG )m 2 n 2 +Q GG (m 4 +n 4 ) 

Q44 * Q4 4m 2 +Q E bH 2 
$4 6 = (Q4‘4“Q6 5 )mn 
$bs * Q4 4 n 2 +Q g 5m 2 

where 

m « cos0 n = s?n0 


and 9 Is the angle between x-axIs and x^axfs measured from 
x to x, counterclockwise as shown In Figure 2.1. 


2.1.2 Plate Strain-Displacement Relations 

The displacement components of the laminated plate are 
assumed to.be of the form [16] 

u(x,y,z) - u°(x,y) + z<£ x (x, y) 

v(x,y,z) = v°(x,y) + z0 v (x, y) (2-6) 

w(x,y,z) = w° (x,y) = w(x,y) 

where u d , v° and w° are the midplane displacement components 
In the x— , y- and z-d l rect i ons , respectively, and <f> x and <f> y 
are rotations of cross-sections perpend I cul ar to x- and y- 
axis, respectively (see Figure 2.2). In Equation (2.6) we 
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Figure 2.2 Laminate displacement components for a cross- 
section perpend I cul ar to the y-axls 
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have assumed that u and v vary linearly In the thickness 
direction, while w Is constant through the thickness. 

The strain components for a point In k-th lamina of the 
laminated plate with a distance z from the midplane can be 
computed as 

® x x k ** ®x° ZKx 
6yy k « 6y° + ZKy 

Txv k - *xy°+ z/c xv (2-7) 

i yx u = 9w/3y + 9v/9z *= 9w/9y + <£ v « Ty* 0 
’V xz k = 3w/9x + 9u/3z =* 3w/9x + «£ x *= Vxz 0 

where 

i x ° = 9u°/9x 

7y 0 « 3v°/9y ( 2 - 8 ) 

T xv °= 9u°/9y + 9v°/3x 

are the In-plane strain components of mldplane,- and 
/c x = 9</> x /9x 

K y = 90 y /3x (2-9) 

*xy = 90x/9y + 9^> y /9x 

are the rotation gradients. 

In Equation (2-7), since w, 0 X and <£ v are Independent of 
z, It follows that the transverse shear strains are constant 
through the thickness of the plate. 
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Equation (2-7) can be written In concise matrix form as 



( 2 - 10 ) 


Thus, the strain components at any point in the plate can be 
determined from the extenslonal strain components of the 
mldptsne, the rotation gradients of the plate and the 
distance z from the ml dpi ana. 


2.1.3 Stress-Resultants and Laminate Constitutive Equations 

Substitution of Equation (2-10) In Equation (2-4) gives 
the stress components for a point In the k-th lamina as: 


k 
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( 2 - 11 ) 


The stress-resul tants acting on a laminate can be 
obtained by Integration of the stresses In each lamina 
through the laminate thickness. Specifically, the In-plane 
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dz 


( 2 - 12 ) 


the stress couples are given by 




ph / 2 
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and the transverse shear forces are given by 


13 ■ i- 


h / 2 
h/ 2 



(2-14) 


The sign convention for these stress-resul tants along 
with the geometry of a typical N-layer laminated plate are 
shown In Figure 2.3. 

Substituting Equation (2—1 1 ) Into the right hand sides of 
the above three equations and performing the ’Integrations, 
we obtain 
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where 

(A, j ,B| j ,D, j) * f h/2 ^i j (1 ,z,z 2 )dz 1,J ■ 1,2, { 

v -h/ 2 

and 


a *u = X. h/2 ^'J dz '-j = 4 - 5 


Equations (2-15) through (2-17) are usually 
symbol leal 3 y as 


'N 


> 

00 

O 

1 


e 0 ' 

M 

• S3 

B D 0 

■ 

K 

.Q. 


1 

* 

< 

0 

0 

1 


i . 


(2-15) 


(2-16) 


(2-17) 


(2-18) 


(2-19) 


written 


( 2 - 20 ) 


which Is the laminate constitutive equation with trantverse 
shear effect Included. 



2.1.4 Plate Equations of Motion 
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The stress-equations of motion for the k-th lamina are 
given by 


XXIX 

+ 

^xyiy 

+ 

Tx z » x 

m 

• * 
P u 

K V'X 


(T y y , y 

-f 

T y z * x 

m 

• t 

pv 

XZIX 

+ 

r V2»V 

+ 

&z z • z 

bn 

pw 


where p Is the mass density, Integrating Equation (2-21) 
through the thickness of laminate and then substituting 
Equation (2-12), (2-14) and (2-^6) In, we obtain 


Nk.k 

+ 

^X¥*V 

n pu° 

+ R^x 


Nx y * x 

+ 

Ny , y 

BS p(/® 

+ R#y 

(2-22) 

Qx » X 

+ 

Qy p y 

+ q ® 

Pw 


where q 

I s 

the 

normal 

traction on the plate. 

Multiplying 


the first two equations of Equation (2-21), Integrating 
through the thickness of laminate and then substituting 
Equations (2-13), (2-14) and (2-5) In, we obtain 

M x ,x + M*v*y “ Qx “ Ru° + I# x 

(2-23) 

M xv ,x + M v , v - Q v * Rv° + I* v 
In which P, R and I are defined as 


(P,R, I) = f h/2 p( 1 ,z,z 2 )dz 
J-h/a 


(2-24) 


Equations (2-*22) and (2-23) are the plate equations of 
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motion. Substitution of Equation (2-20) and then the 
strain-displacement relations In these two equations yield 
the equations of motion In terms of midplane displacements 
and rotations of the plate, 

A graphite/epoxy laminated plate provided by NASA Lewis 
Research Center was used throughout this study. This 
laminate Is a [0°/45 0 /0 0 /-45 0 /0° J 2 s graph! te/epoxy composite 
with 0.0053 Inch ply thickness and the following ply 
properties El 5]: 

E 1 = 17.5 X 10 6 psl. 

E 2 = 1.15 X 10® psl . 

0 1 2 25 Gis = G 23 = 0.8 X 10 6 psl. (2-25) 

v i 2 — 0 . 30 

p = 0.000148 1 b-sec 2 / I n 4 

For symmetr leal 1 y laminated composite plate, B fJ =0 and 
R =0. In addition, by choosing the x-axIs of the laminate 
reference system to coincide with the 0 ° fiber direction, we 
obtain A t 6 = A 26 - 0 and D, 6 = D 26 . Further, In this study, 
we assume G 13 - G 23 = G 12 , and consequently, A* 4B = 0 and 
A* 4 4 = A*s s • For this particular laminate, the 

displacement-equations of motion are given by 

A, i 9 2 u°/9x 2 + A G 6 9 2 u°/3y 2 + (A , 2 + A 6 6 ) 9 2 v 0 /9x9y = Pu° 


(A , 2 + A 66 )9 2 u°/9x9y + A 66 9 2 v°/9x 2 + A 22 9 2 v°/9y 2 = Pv° 
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D n 9 2 ^/ax 2 + 2D ie 9 2 0 x /9x9y + D 66 9 2 <*> x /0y 2 

+ 0, 6 (9 2 0 y /9x 2 + 9 2 0 y /9y 2 ) + (D 1a + D 6 «)9 2 <A v /9x9y 
-A* 44 (9w/9x + 4 > h ) » I* K (2-26) 

D 16 (9 2 ^ x /9x 2 + 9<^ x /9y 2 ) + (D ia + D 6 6 ) 9 2< ^x/9x9y 
+ D 66 9 2 0 v /9x 2 + 2D 16 9 a ^ v /9x9y + D a a 9 2 </> y /0y 2 
-A* 44 Ow/9y + 0y) “ I# v 

A* 44 (9 2 w/9x 2 + 9 2 w/9y a + 9<£ x /9x + 9^ v /9y) + q ** Pw 

In Equation (2-26), the first two equations govern the 
In-plane motion while the last three equations govern the 
flexural motion. 


2.2 Propagation of Harmonic Waves 

Consider a Infinitely large laminated plate governed by 
the equations of motion (2-26). We assume plane harmonic 
waves In the form 

u° = U exp[ I k (r? - ct) ] 
v° ® V exp[ I k ( 77 - ct) ] 

w = W exp [ 1 k ( 7? - ct)] (2-27) 

0* = $ x exp [ 1 k ( 7? - ct)3 
<£ v « $ v exp[lk(7? - ct) ] 

propagating over the plate, where U, V, W, 3> x and iS y are 
constant amplitudes, k Is the wave number, c Is the phase 
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velocity and >7 Is given by 

V m x cqsk + y slnoi 

In which oi Is the angle between the direction of wave 
propagation and x-axls, 

Substitution of Equation (2-27) into Equation (2-26) wtth 
q *» 0 yields a system of five homogeneous equations for the 
five constant amplitudes. In order to have a nontrivial 
solution, tho determinant of the coefficient matrix Is set 
equal to zero, Since the equations are uncoupled Into two 
groups, the determinants! equation can be separated Into two 
equations os 

|fi u | w 0 (2-29) 

for the In-plan© extenslonal and In-plan© shear waves, and 

jb, j | * 0 (2-30) 

for the flexural waves. In Equations (2-29) and (2-30) the 
coefficients a u and are given by 

* A n cos a o! + A Bfl sln a a - Pc a 

fi, j « a a t * (A^a + A& # ) s ! ruxcostt (2-31) 

a aa ** A 6fi cos a o; f A aa sln a ot - Pc a 

and 

b n » D 1i k a cos a o: + 2D, dk a slna:cosa + D e6 k a sln a o{ 

+ A * 44 Ik a c a 


ofpqq *}■ base ® 

h °or quality 


(2-28) 
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b, a * bai » D 1 f k a co« J? o( + (Dta + D* t )k a s I nacosot 
+ D,ik a $ln a o< 

bia « b 81 * I A* 4 4 kcoso( (2-32) 

baa - D ae k a cos a o! + 2 D, <$k a s t nwcoscx + D aa k a sln a « 

4* A% 4 - Ik a c a 

bja m b 9 a » i A% 4 ks I na 

b 35 « “A* 44 k a 4 - Pk a c a 

Expending Equation (2-29) wo obtain a quadratic equation 
In c a as 

c 4 - d,c a 4* d a * 0 (2-33) 

where 


d, * (A 1t cos a oi 4* Aaasln a « + A 8 b)/P 


da 


A n cos a <x 4* A« 6 sln a a (A, a * A* c ) s l notcosot 
(A, a * A e e ) s I noicosa A 66 cos a a 4- A aa sln a ci! 


(2-34) 


It Is noted that the phase velocity c does not depend on 
the wave number k, thus these waves are nondlsperslve. In 
studying of transverse Impact problem where In-plane 
deformation Is negligible, this nondlsperslve property has 
no significant effect . Should In-plane deformation become 
Important, higher order approximation of displacement 
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components u and v must be assumed and the dispersive 
property of these In-plane waves could be Included. 

From Equation (2-34) It Is evident that there exist two 
phase velocities correspond I ng to two modes of wave. 
Although these two waves Involve both In-plane extenslonal 
deformation as well as In-plane shear, from the eigenvectors 
we are able to tell which one Is dominant. Thus we label 
the two waves as fn-plane extenslonal wave and In-plane 
shear wave accordingly. 

The determl nantal equation given by Equation (2-30) 
yields three positive roots In c 2 Indicating that three 
flexural waves exist. These phase velocities are functions 
of the wave number k, thus they are dispersive. Among these 
three modes of wave, only the lowest one correspond! ng to 
the transverse shear wave has a finite ^velocity as k-*0 or as 
the wave length becomes Infinite. The dispersion curves for 
the waves of the lowest mode propagating In the directions 
of 0°, 45° and 90° respectively are plotted In Figure 2.4 
with the non-d I mens lonal phase velocity vs. the non- 
dimensional wavelength A/h. It can be seen that they all 
approach the value of VG 13 /p as the wavelength becomes 
shorter. The phase velocities for the two higher modes, 
however, approach different values In different propagation 
directions when \-0. For laminated composite which are 
anisotropic In general, the phase velocity varies from one 
direction to another. As a result the wave surface will 
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Figure 2.4 Dispersion curves for plane harmonic waves 
propagating In the O 0 - 45°- and 90°- 
dlrect I ons 
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become a rather complicated shape as it propagates. This 
will be discussed In the next section. 

Substitution of w b kc In Equation (2-32) yields a set of 
frequency equations for flexural waves. Figure 2.5 shows 
the frequency curves of these waves for a « 0° , 45° and 90°, 
respectively, with the non-dimensional frequency vs. the 
non-dimensional wavelength. The cutoff frequencies for the 
two higher modes have a value of yi2Gi 3 /p/h. Comparing with 
the exact cutoff frequency < 7t/h ) \/G t 3 /p lt It can be seen that 
If the shear correction factor rr 2 / 1 2 Is Introduced, this 
theory will predict the correct cutoff frequency. 


2.3 Propagation of Wave Front 

Impact of foreign objects on a laminated plate with a 
very short duration could generate weak shock waves which 
will propagate Into the rest of the structure with finite 
velocities, and the positions of the wave fronts define the 
regions being disturbed at any Instant after Impact. 
Damages to the laminated plate may possibly occur as the 
first wave front hits the weakest part. It Is hence 
Important to Investigate the propagation of these shocks In 
the plate. There have been works dealing with the 
propagation of wave front In anisotropic elastic media [20- 
22]. Moon [23] presented an analysis of wave surfaces In a 
laminate by treating It as an equivalent homogeneous 
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orthotropic plate. The acceleration waves and their wave 
fronts were Investigated. The propagation of shock waves In 
more general laminates which exhibit the bendlng-extenslonal 
coupling were studied by Sun [2] * The ray theory was 
employed to construct the wave front surface. The growth 
and decay of the shock strength were also discussed. In 
this section, the analytical procedures developed by Sun [23 
were applied on a [0°/45 0 /0 0 /~45 0 /0°3 as graphlte/epoxy 
laminated plate. 


2.3.1, Kinematic Conditions of Compatibility on the Wave 
Front 

A wave front, which will be denoted by n, fs defined as a 
surface travelling over the plate as time varies 
continuously, and across which there may exist a 
discontinuity In the stress, particle velocity and their 
der I vat I ves . 

Consider a discontinuous surface Q passing some 
observation point In a medium at a certain time t, Let f“ 
be the value of a field function f(x,,t) (e.g. stress, 
particle velocity, etc.) behind the surface 0, and f + be 
the value of f In front of It, then the discontinuity of 
function f can be expressed as 

[ f j « f+ - f- (2-35) 
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In which the right hand side Is to be evaluated at the time 
and location on n passing the observation point, and the 
Jump across the wave front Is denoted by square bracket. 

Surface n may be expressed mathematically by an equation 
of the form 


sli(x, ,t) « 0 (2-36) 

or, by making t explicit, as 

ili(x ( ,t) - F(X| ) - t - 0 (2-37) 

which represents a family of surfaces In x, -space with t as 
a parameter. By evaluating f + and f” at t » F(xi), the Jump 
of f across the wave front becomes 

Cf (x, )] - f + (x, ,F<x, ) ) - f" (x, »F(x, ) ) (2-38) 

The rate of change of Cf3 for an observer moving with ft 
Is given by 

d[fj/dt » (8f + /0x, - ef-/dxi )dx,/dt + (9f + /et - df-/dt) 

« c, [9f/0X| ] + E0f/0t3 (2-39) 

where t = F(x, ) Is substituted, and c ( <= dX|/dt are velocity 
components of wave front relative to the material. 

If the laminate theory Introduced In previous section Is 
used, then the plate displacement components are u° , v°, w, 
and <£ v , while the spatial variables are x, = x and x 2 - 
y. It Is assumed that the Integrity of the material Is not 
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affected by the propagation of the stress wave front, then 
these displacement components will remain continuous. 
Consequently, we have 

[u°] * i v°3 * M * M - l<t > V 1 - o (2-40) 

across the wave front. Applying the- general condition of 
Equation (2-39) on u° , together with Equation (2-40), we 
obtain 

[au°/3xj3cj + [u°3 « 0 J * 1,1 (2-41) 

Let c n and nj be the normal velocity and the unit normal 
on the wave front, respectively, It follows that 

n j c j «* c n (2-42) 

and Equation (2-41) becomes 

Cau°/3xj3 - -Cu°] nj/c n J -1,2 (2-43) 

Similar relations can be derived for the other 
displacement components v°, w, and <£ v . Together they 
specify the kinematic conditions of compatibility on the 
wave front. 


2.3.2 Dynamical Conditions on the Wave Front 

Consider a finite volume V of a material medium and 
denoted by S the boundary or surface of V. For a continuous 
and differentiable function f(X|,t) In V, It can be shown 


[23] that 
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^ v f(x, ,t)dV - J v f,tdV + J # GfdS (2-44) 

under deformation of the medium, where G Is the normal 
velocity of the surface S, In the case where the 
deformation of the volume V Is determined solely by the 
motion of material particles, we have 

G * u,n ( * u n (2-45) 


where u ( Is the displacement components, n, Is the outward 


onrl 

Wl IM 


W1 tv* 


nAr*mt!s I 

» |WI HIM » 


1 a#*» t fu 

t v» i vv* i w y 


of 


particle on S. If there exists a discontinuity surface (or 
wave front) travelling with velocity c, In the medium, by 
choosing this surface as the boundary of V, we have 


G « c | n | « c n 


(2-46) 


where c n Is the normal velocity of wave front. 


Suppose that a volume V whose motion Is determined by the 
deformation of the material medium, Is divided by a 
travelling surface Q Into two volumes V" and V + as shown In 
Figure 2.6. The surface S Is also divided Into two portions 
S“ and S + which form parts of the boundaries of V" and V + , 
respectively. The remaining part of the boundary Is formed 
by n 0 which Is a segment of n. In Figure 2.6, n, denotes 
the unit normal of fl In the direction of travelling, and ni* 
denotes the unit outward normal of S. 
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Taking f " p u t In Equation (2-44) and using aquation (2- 
45) and (2-46), we obtain 


4f v J> Cl T dV - T>.tdV + J’ f U„-pi 1 7dS + J^CnpOTdO (2-47) 

^puldv - JjpuD.tdV + fj>*pV,<tS - J^o^utdn (2-40) 


where u~ and u| are the velocity components of material 
particle In V" and V*, respectively, Combining the above 
two equations gives 

alJ/MV - J v (pO,), t dV + J^pujdS + J^CiJpC J|d$ 

+ f c n p(uT - 0|)dO (2-49) 

v n* 

From theory of elasticity w® have 


^/ v PU,dV - J o cr,jnjdS (2-50) 

If we let the volume V approach zero at a fixed time In such 
a way that It will pass Into fl 0 * then the volume Integral In 
Equation (2-49) will evidently approach zero; however 


JViJpuldS - - /jSJpiifdn 


(2-SI) 


J^Ci-pCi-dS - Jj^pCiTdQ 
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(2-52) 


J o cr U RjdS ~ 


^<'^4 “ tfTj)njdn 


(2-53) 


whore <r7j and cr| j are the stress components on the sides of 
O 0 , respectively. 

Substituting Equations (2-50) through (2-53) Into 
Equation (2-49) gives 

P (tf| .-cry j )njdfl = f jouJ (c n -u“ )dfi - P pu| (c n -u?;)dn (2-54) 

tlno J J t/fto’ v nt r 

Using Co* j j ] and CCi|] to represent the Jumps of stress and 
particle velocity across the wave front, and utilizing the 
fact that c n » u n , we obtain 

f tff,j3njdn « - f pc n Cu,3dn (2-55) 

v ft® v no 


Since this condition Is Independent of the extent of the 
surface Integration Q 0 , It follows that 

Ccr | j 3 n j « - pcnCCi,] (2-56) 

In the case of laminated plate, 1 = x,y,z and J « x,y. 

Substitution of Equation (2-6) Into Equation (2-56) 


yields 
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Ctrl j 3 nj - « pc n { Cu°] + zt<M> 
[tfajlnj - - pc n { Cv°] + zE<M> 

C^ajlnj « - pc n [w] 
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Integrating Equation (2-57) over the thickness of plate 
gives 


[N x 3n x + EN xy 3n y = - Pc n [u°3 - Rc n [<fc x ] 

CN xy 3n x + [N v ]n y = - Pc n Ev°] - Rc n [£ y ] (2-58) 

CQ x 3n x + [Q y 3n y = - Pc n [W] 

Multiplying the first two equations of Equation (2-57) by z 
and then integrating over the thickness, we obtain two more 
equat I ons 

[M x ]n x + [M xy ]n y = - Rc n [u°] - Ic n [£J 

(2-59) 

[M xy ]n x + [M y 3n y = - Rc n [v°] - Ic n [£ y ] 


where P, R and I have been defined In Equation (2-24) 

The five equations given by Equations (2-58) and (2-59) 
are the dynamical conditions on the wave front for the 
laminated plate, 


2.3.3 Propagation Velocity of the Wave Front 

Across the wave front, the laminate constitutive 
relations given by Equation (2-20) can be written as 
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(2-61) 

are the Jumps of the stress resultants, and 

{ [e] } T = { £du°/ax3 , E9v°/9y] , Cau°/dy] + [av 0 /ax] } 

{ C/cl } T « {[a^/axi.^v/ayi.Cd^/ayj+ce^/dx]} (2-62) 

iM) T - { caw/9y] , [aw/3x3 } 

are the jumps of the strain components. In Equation (2-62), 
the conditions [0 X 3 = [<£ y ] « 0 are substituted. 

Substituting of Equation (2-43) and the similar relations 
for other kinematic variables In Equation (2-60), we can 
express the jumps of the stress resultants In terms of the 
Jumps of the time derivatives of the displacement components 
u°, v°» w, and 4> v . These relations are then substituted 
In Equations (2-58) and (2-59), which results In five 
homogeneous equations. For [0°/45°/0 0 /-45 0 /0 0 ] 2S graphite/ 
epoxy laminated plate which Is symmetrical and balanced 
(I.e. B | j — 0 , A j 6 = A 2 B = 0 , R = 0 and D ^ g — D 2 $ ) » these 
five equations are uncoupled Into three groups as 


CN3 
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CD 
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CM3 

• S3 

B D 0 
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U3 

EQ3. 
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where 

{[N]F - { [N x 3 , [N y ] , [N xv ] } 
{ EM] } T - { [M x 3 , CM V 3 , [M x v 3 } 
{ CQ3 } T - { E r *3 , CQ V ] } 
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Ca.j] 

tb tJ ] 



M 0 

m 0 


(A* 4 4 - PGn a )C\V 3 » 0 


(2-63) 

(2-64) 

(2-65) 


In which Ca t j 3 and [b|j] are both two by two symmetric 
matrices, and their entries are given by 


a, , 

m 

n* a A n 

+ 

n Y 

a A*o - Pc„ a 





a 

m 

«ai ■ 

n«n 

i v (A\ a 4* Afce ) 



(2-66) 

®aa 

m 

n* a A« 6 

4* 

n Y 

a Aaa - Pc„ a 





bi i 

m 

n x a D n 

+ 

2n 

xn v D 1a 4* n v 

a D«# ** 

ICn 3 



bi a 

m 

ba i * 

Di„ 

4* 

n x n v (D ia 4- 

Dq q ) 


(2-67) 

ba a 

m 

n x a D«« 

4- 

2n 

*n v D 1a + n v 

a ^aa ~ 

ICn* 




It can be seen that Equation (2-63) describes the In- 
plane extenslonal and the ln-plane shear wave fronts, 
Equation (2-64) describes the bending moment and the 
twisting moment wave fronts and Equation (2-65) describes 
the transverse shear wave front. 

From Equation (2-65), we obtain the normal velocity with 
which the transverse shear wave front propagates as 

c„ a - A*4 4 /P (2-68) 

It Is noted that this velocity Is Independent of the 
direction of propagation, and Is called directionally 
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Isotropic wave front. 

Equations (2-63) and (2-64) yield non-trtvlal solutions 
only If the determinant of the coefficients matrices vanish, 
l ,e, 


| a i j | "> 0 (2-69) 
| b, j | - 0 (2-70) 

Each of the above equations can be expanded Into a 
quadratic equation of c n a , For [0°/45°/O a /-4S 0 /O 0 ] as 
graph l te/epoxy laminated plate, the normal velocities of 
wave fronts corresponding to the ln-plane modes and flexural 
modes are plotted In Figure 2,7 and 2.8, respectively. It 
Is noted that the normal velocities of the ln-plane 
extensions! and ln-plane shear modes are symmetrical about 
x-axls and y-axls, while there Is no such symmetry for the 
bending moment and twisting moment modes, 


2,3,4 Wave Surface and Ray 

From Figure 2.7 and 2.8, It can be seen that for 
laminated composites which are anisotropic In general, the 
ln-plane and flexural wave fronts travel with different 
normal velocities In different directions. In other words, 
the Initial shape of a wave surface will be distorted after 
It propagates, However, Equations (2-66) and (2-67) show 
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Figure 2.8 Normal velocities of flexural wave fronts 
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that for any fixed normal direction nj , c n Is a constant. 
Connecting the points having the same unit normals to the 
travelling wave front surface, we obtain a family of lines 
which are called rays. Thus, along a ray, the normal 
velocity of wave front remains unchanged. By using the ray 
theory which has been well established in the field of 
geometrical optics, we are able to construct the wave front 
surface , 

Recall Equation (2.37) 

F(X| ) - t = 0 1*1,2 (2-37) 

which represents a family of wave fronts propagating over 
the plate with t as a parameter. It follows that 


dF/dt « OF/3x,)(dx,/dt) - (3F/3x,)c, * 1 (2-71) 

By putting 

p, = 3F/3x t * VF (2-72) 

Equation (2-71) becomes 

p , c | * 1 (2-73) 

Since p, is normal to the surface F, It can be written as 
p, = | pi] n, (2-74) 


where jp, | denotes the length of p ( . Combining (2-73) and 
(2-74), we obtain 
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i P i I n | c , « | p » | o n * 1 (2-75) 

from which we obtain 

p, - n,/c n (2-76) 

In Equation (2-76), pi Is called the slowness vector 
which has the direction normal to the wave front with the 
magnitude being equal to the Inverse of normal velocity c n . 

Subs 1 1 1 tutlng Equation (2-76) In Equation (2-69) and (2- 
70>, we obtain two equations In terms of p, 

P* 2 An + p ¥ 2 Ass “ P PxP v (A 13 + A 66 ) 

m Q 

PxPyCAia + A 6 e) Px 2 Aee + Py 2 Aaa “ P 
P x 2 U i i +2p x Py D ^ s+Py a 0c g — I D 1 6 +p x Py ( D i a+D G e ) 

* 0 

Di 6+PxPy <D t a+D oe ) p x a D 66 +2p x p v D 1 G +Py 2 Daa“I 

which can be written In a general form as 

g(p, ) “0 1*1,2 (2-77) 

In view of Equation (2-72), we recognize that Equation 
(2-77) may be regarded as a set of first-order partial 
differential equation for F. A standard method of solving 
first-order partial differential equation Is by means of 
characteristics [24], which reduces the equation to a system 
of first-order ordinary differential equations. In our 
case, Equation (2-77) then Is equivalent to the following 
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dx/ds * 9g/9p x dy/ds * 9g/9p y (2-78) 

dp x /ds * -9g/9x dp y /ds * -9g/9y (2-79) 

where s Is a parameter, These equations, together with 
Equation (2-77) describe the ray geometry and the normal 
dt recti on of the wave front propagating along the ray. 

From Equation (2-78), we have 

dy/dx “ ( 9g/9p y ) / (9g/9p x ) (2-80) 

Since the normal direction of wave front along a ray Is 
constant, It can be seen from Equation (2-76) that p, Is 
also constant along a ray. For laminated composite which Is 
assumed to have homogeneous material properties, Equation 
(2-77) shows that g(p,) does not depend on x, , consequently, 
9g/9p x and 9g/9p y are all constants along a ray. Thus, the 
solution of Equation (2-80) Is then given by 

y « f(x - x 0 ) + y 0 (2-81) 

where x 0 and y 0 are the Initial values of x and y at t = 0, 
and f « ( 9g/9p v ) /( 9g/9p x ) . This equation shows that the 
rays In a homogeneous solid are straight lines. 

From Equations (2-73) and (2-77), we have 

C| dpj w 0 (2-82) 


dg « (9g/9p, ) dp, « 0 


(2-83) 
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Eliminating dpi from those equations yields 

dx,/dt « c, « Og/ap, )/(pj9g/apj) (2-84) 

where summation over J Is understood. 

Equation (2-84) can be solved to obtain the position of 
wave front at time t. Again, since 9g/9pi and pi are all 
constant along a ray, we obtain the solution of Equation (2- 
84) as 

X = (ag/9p x )t/(pj9g/9pj ) + x 0 (2-85) 

Y — Vvvj/Vfjy/Wf vpj vy/ W(^j / * Y O \ vU/ 

where x 0 and y 0 denote the Initial wave position at t ■ 0. 
Consider at t = 0, a wave front forms a circle given by 


x 0 = h cosod 
y 0 = h si nod 


(2-87) 


At this Instant, the normal directions to the wave front 
coincide with the radial directions. Due to the different 
velocities of propagation in directions, this Initial shape 
would be distorted. By using Equations (2-85) and (2-86), 
the subsequent positions of the wave front can be 
determined. Figures 2.9-2.12 show the wave front positions 
at two consecutive Instants after t = 0 for the In-plane 
extensional, In-plane shear, bending moment and twisting 
moment modes, respectively, for the [0°/45 0 /0 0 /-45°/0 0 ] 2S 


ORIGINAL PAGE \S 
OF POOR QUALITY 


42 


graph 1 to/©poxy laminated plate. It Is noted that for 
symmetrical laminates, the In-plane modes are uncoupled from 
the bending modes. The rays along which the normal 
directions to the wave front are 0°, 45° and 90°, 
respectively, are also shown In the figures. It Is seen 
that the wave fronts of the In-plane extenslonal and the 1 n~ 
plane shear modes possess symmetry With respect to x-axIs 
and y~axls. The wave fronts of the bending and twisting 
moments, however, lose their original symmetry with respect 
to x~axls and y-axls, This Is an Indication that In 
performing analysis of flexural deformation of this 
laminate, one can not take a quadrant for analysis, a 
practice followed by many authors dealing with homogeneous 
and Isotropic plates. 

From Figures 2.9-2.12, It Is also Interesting to note 
that ray geometries for these two groups of wave fronts are 
quite different. For the In-plane extenslonal and In-plane 
shear wave fronts, the rays coincide with the normal 
directions when a ® 0° and 90°. Along other directions, the 
direction of the ray deviates from the normal direction of 
the wave front. It was discussed In [2] that the degree of 
spreading of rays Is proportional to the decay of the stress 
amplitude at the wave front. Thus, from Figures 2.9 and 
2.11, one can conclude that the strength of the in-plane 
extenslonal and bending moment wave fronts decay more 
rapidly In the y-dlrectlon than In the x-dl recti on. 
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A photoelastfc study of anisotropic waves in a fiber 
reinforced composite has been done by Dally et aK [9], 
The waves was produced by a explosive charge In a small hole 
on the plate. The result showed clearly an elliptic-like 
stress wave front pattern. This indicates that stress waves 
In anisotropic materials propagate with different velocities 
In different directions. 
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CHAPTER 3 

STATICAL INDENTATION LAWS 

A brief Introduction of the historical development on 
Impact problem Involving homogeneous Isotropic materials was 
given by Goldsmith [123. Hertz [11] was the first to obtain 
a satisfactory solution on contact law for two Isotropic 
elastic spherical bodies. When letting the radius of one of 
the spheres go to Infinity, this law then describes the 
contact behavior between a sphere and an elastic half-space. 
The Hertzian law, In spite of being static and elastic In 
nature, has been widely applied to Impact analyses where 
permanent deformations were produced. The use of this law 
beyond the elastic limit has been Justified on the basis 
that It appears to predict accurately most of the Impact 
parameters that can be experimentally verified. 

In studying Impact responses of laminated composites, the 
problem becomes extremely complicated. One may easily 
realize that the Hertzian contact law which was derived 
based on homogeneous Isotropic materials may not be adequate 
In describing the contact behavior of laminated composites 
due to their anisotropic and nonhomogeneous properties. 
Moreover, most of the laminated composites have finite 
thickness which can not be represented by a half-space. In 
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many existing analytical works [253, loadings to the 
1 amt nates w era assumed known, and the responses of the 
laminates were assumed elastic. 


Will Is [26] obtained expllct' formulas for Hertzian 
contact law for transversely Isotropic half-space pressed by 
a rigid sphere, and extended It to the application of Impact 
problems. It was shown that 

F « ka" (3-1 ) 


with n « 3/2 Is valid for the contact force F and the 
Indentation ci, where k Is a contact coefficient whose value 
depends on the material properties of the target and the 
sphere, and the radius of sphere. 

A modified contact law with 


R 1 f 2 

k = (4/3) - S (3-2) 

1 - y,s 2 

E s E t 


was used [13] In an analytical study on Impact of laminated 
composites. In Equation (3-2), R s ,v s and E s are the radius, 
Poisson's ratio and Young's modulus of the sphere, 
respectively, and E t Is the Young's modulus of the laminates 
In thickness direction. It was also suggested by Sun et. aj_. 
[27] that the value of k can be experimentally determined. 
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Recently Yang and Sun [14] have conducted static 

Indentation tests on the [O 0 /45 0 /O 0 /-45°/O 0 ] 2S graphite/ 

epoxy laminates using spherical steel Indenters of 0.25 In. 

and 0.5 In. diameters. The results were fitted Into 

Equation (3-1) and were found that the 3/2 power Is valid. 

In addition, It was also observed that even for small 

amounts of load there were significant permanent 

Indentations.. This Implies that the unloading curve has to 

be different from the loading curves.. In order to account 

for the permanent deformation, the equation 

/ot - ec e \ q 

F 3 F m \ (3-3) 

\a m ~ a 0 / 


proposed by Crook [28] was used to model the unloading path 
where F m Is the contact force at which unloading begins, a m 
is the Indentation corresponding to F m , and a 0 denotes the 
permanent Indentation In an unloading cycle. Equation (3-3) 
can be rewritten as 

F = s(a - a 0 ) q (3-4) 


in which 


s = F m /(a m - a 0 ) q (3-5) 

Is called unloading rigidity. In order to simplify the 
modeling of the unloading law, It was assumed [14] that the 
value of s for all the unloading curves remains the same. 
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Consequently, a constant a cr given by 

a cr * k/s (3-6) 

was Introduced. It was also shown that q»5/2 fitted the 
unloading path very well, and the permanent Indentation a 0 
was then related to a m by 

a 0 /«m 82 1 - (od cr /a m ) a/B as > ct cr 

(3-7) 

a 0 - 0 as £ a cr 

The value of a cr was found to be Independent of the size 
of the Indsnter and hence can be regarded as a material 
constant . 

It was also mentioned In [14] and [29] that there were 
some practical difficulties In performing the tests. Since 
the Indentation was measured step by step using a dial gage 
and readings on the gage were taken about 10 to 20 seconds 
after the load was Increased by one step, the creep effect 
may cause an appreciable error to the results. Another 
Important problem was that It was almost Impossible to 
measure the permanent Indentation accurately using the dial 
gage. In order to overcome these problems, a Linear 
Variable Differential Transformer (LVDT) was used In this 
study to measure the Indentation. 

The LVDT Is an electromechanical transducer that produces 
an electrical output proportional to the displacement. 



52 


ORIGINAL. PAGE ES 
OF POOR QUALITY 

Connecting this output and the one frcm the strain Indicator 
which Is used to measure the applied loading to a X-Y 
!p Jotter, one can obtain a continuous loading-unloading 
curve. By changing the loading rate which can be applied as 
fast as 50 lb. /sec., It Is possible to examine the 
Significance of creep effect on the contact law. The 
starting point and final point of a loading-unloading cycle, 
which represent respectively the Instants of contact and 
separation of the Indenter and the specimen, can be easily 
determined from the curve. Thus, the measurements of 
permanent Indentations are much more accurate than those 
using the dial gage. 


3.1 Specimens and Experimental Procedure 

Two groups of test specimens were prepared from a [0°/ 
45°/0 0 /-45°/0 0 ] as graph l te/epoxy laminate. They were cut In 
the way such that the longitudinal axis of the beam specimen 
of the first group was parallel to the 0° fiber direction 
while the second one was perpend Icul ar to It. The latter 
then becomes [9O°/45 o /9O 0 /-45 o /9O° 3 as laminated beams. The 
thickness of the beam was 0.106 In. and the width was 
approximately 1.25 in.. In all tests, the specimens were 
clamped at both ends. It was shown In [14] that the span of 
the specimen In the range of 2 In. to 6 In. has little 
effect on the contact law. Hence, only one span, l.e. 2 
In., was used In the test, 



ORIGINAL PAGE IS 
OP POOR QUALITY 


53 


The experimental set-up Is shown schematically In Figure 
3.1, LVDT was mounted on a 'C' bracket fixed to the loading 
piston so that only the relative movement between the 
Indenter and the specimen was recorded. The load was 
applied pneumatlcal 1 t by a plunger and It was measured using 
a load cell and a strain Indicator. Outputs from LVDT and 
strain Indicator were fed Into an X-Y plotter so that a 
continuous force- Indentation curve can be obtatned. Two 
spherical steel Indenters of diameters 0.5 In. and 0.75 In. 
were used. 


3.2 Experimental Results 
3.2.1 Loading Curves 

The experimental curves were first digitized Into some 
discrete data points and then fitted Into Equation (3-1) 
using least-squares method. Figures 3.2 and 3.3 show the 
test data and the fitted curves for 0.5 In, diameter 
Indenter. It can be seen from these figures that the 3/2 
power Index gives very good results. However 5 the contact 
coefficient k of [0°/45°/0 0 /-45 o /0 0 ] 2S specimen Is less than 
the one of [90 0 /45 0 /90 0 /-45°/90 0 3 as specimen by about 7 %. 
During the test, larger deflections were observed for the 
second group of specimen due to their lower flexural 
rigidity. This means that the contact area Is also larger 
and the Indentation under same amount of loading should be 










K=1 .2838E+06 





K=1 .3756E+06 
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smaller comparing with the first group of specimens. 
Consequently, the higher value of k for the [90°/45°/90 0 /- 
45°/90°] 2C specimens Is reasonable, 

The results for 0.75 In. diameter Indenter ore presented 
In Figures 3,4 and 3,5. Again, good agreement between the 
experimental data and fitted curves Indicates that the 3/2 
power Index for loading law Is valid. The values of k for 
both Indent^rs are summarized In Table 3.1, It should be 
noted that the average value of k obtained from the two 
groups of specimens was used later In a finite element 
analysts of Impact responses. 

3.2.2 Unloading Curves 

By choosing a suitable value for q, It can be seen from 
Equation (3-5) that once the relation between c t 0 and oi m Is 
established, the unloading rigidity s Is then determined. 
T*‘©t results show that the permanent Indentations a 0 and the 
corresponding maximum Indentations a m exhibit a rather 
linear relationship. The equation given by 

«:i 0 = s p (a m - a p ) (3-8) 

Is obtained from the test data for both 0.5 In. and 0.75 
In. Indenters using least-squares fitting method, and are 
plotted In Figure 3.6. In Equation (3-8), ot p can be 
considered as a critical value of Indentation, Once the 
amount of Indentation exceeds u p , permanent deformation will 


occur . 


K=1 ,8326Et06 



Figure 3.4 Load 
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Table 3. 1 

Contact coefficient k of loading law F = ka 1 - 5 


Size of 
Indenter ( I n) 

0.5 

0.75 

Specimen 

Group 1 + 

Group 2$ 

Group 1 + 

Group 2t 

kdb/in 1 * 5 ) 

1 .284x10 s 

1 ,376x10 s 

1 .833x10 s 

1 .990x10 s 

Average k 

1 .330x10 s 

1 .912x10 s 

Ref . [14] 

9.694x10* 



+ [0°/ 45 0 /O 0 /“45°/0 0 ] 2 s specimens 
t [90°/ 45 0 /90 0 /-45 0 /90 0 ] a s specimens 




OOOh 



Figure c 
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Substitution of Equation (3-8) and (3-1) Into Equation 
(3-5) yields 

ka m a/ 2 

s » — — If oim £ a p (3-9) 

CO - Sp)d m + Sp«pj q 

k<x m ? ' 2 

s M If < <x p (3-10) 

^m q 

These two equations along with Equation (3-4) are then used 
to fit the experimental unloading curves In finding the 
value of q. 

Yang [14] has shown that, q * 2.5 fits the test results 
for both 0.25 In. and 0,5 In. Indenters quite well. In 
this study, however, the values of 2,2 and 1.8 were found to 
give the best fitting for 0.5 In. and 0,75 In. Indenters, 
respectively using the aforementioned method (Figures 3,7- 
3,10), For convenience, q = 2,5 was used for 0.5 In. 
Indenter while q = 2.0 was chosen for 3/4 In. Indenter. 
The results, of the curve-fitting are presented In Figures 
3.11-3.14. Further discussions on the unloading law will be 
given In Section 3,3. 
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Figure 3.11 Unloading curves or L0 u /4b u /0'V-4b'VU~j 2 s 
sDeeimens with 0.5 Inch indenter (q=2.5) 
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3,2,3 Reloading Curves 
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The equation 

F » ki (a - « 0 ) p (3-11) 

suggested by Yang [14] was used to model the reloading 
curve, where ki Is called reloading rigidity and p ■« 3/2 was 
found to fit the experimental data quite well. It was also 
observed that the reloading curve always returns to where 
the unloading began, and hence the reloading rigidity can be 
determined by 

ki - F m /(cd m - #«,)*'» (3-12) 

In other words, the reloading test Is not necessary provided 
the unloading condition Is specified. Some reloading curves 
obtained following Equations (3-11) and (3-12), and the 
experimental data are presented In Figures 3,15-3.18, 


3,3 Discussion 

As mentioned before, due to creep the loading rate may 
affect the contact law (l.e. the value of k) , A series of 
tests with different loading rates was performed to examine 
this point. The maximum loading rate the test equipment can 
apply without exceeding it's capacity is about 50 lb/sec.. 
It was found that in the range of 5 Ib/sec. to 50 lb/sec., 
the values of k showed very little scatter, and the effect 
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Figure 3.16 Reloading curve of [90°/45 0 /90 o /-45°/90 o ] 2S 
specimens with 0.E5 Inch Indenter (p=l .5) 
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due to local material nonhomogeneity In the composite may be 
even greater than the one due to the loading rate. However, 
an appreciable decrease of the value k was observed when the 
loading rate was lower than 1 lb/sec.. In some extreme 
cases where loadings were applied as slow as 10 lb/m In., the 
average value of k for 0.5 In. I ndenter was very close to 
the one obtained previously by Yang [14] using dial gage to 
measure the 'Indentation. In this study, the loading rates 
for all tests were approximately equal to 10 lb/sec.. 

Unlike the exponent n of the loading law for which the 
value of 3/2 seems to yield good agreement with all 
experimental data, the exponent q of the unloading law 
(Equation 3-3 or 3-4) reveals much wider deviation for 
different sizes of l ndenter. Value of q = 3/2 correspond! ng 
to an elastic recovery according to the Hertzian theory was 
previously used by Crook [28] In a study of Impacts between 
metal bodies. The experimental results from [14] and 
present study show that the value of q varies from 1.5 to 
2.5. Local plastic deformation, anisotropic properties of 
composite material and unloading rate are all possible 
causes for this deviation. Obviously, an analytical study 
to determine the value of q as function of aforementioned 
factors Is Impracticable. Since the purpose of this study 
Is to establish a contact law that can be used In the 
analysis of Impact, the validity of this law must be 
verified from Impact experiment. This will be Investigated 
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From Equation (3-3) or (3-4), It can be seen that a 0 
plays an essential role In the unloading law and hence the 
value of It must be estimated accurately. Both of Equation 
(3-7) used by Yang [14] and Equation (3-8) used In this 
study for calculating a 0 were obtained experimentally, in 
which a cr and a p are considered to be material constants and 
were determined using a 0 and a m from test data. However, It 

was pointed out In r 1 4] that the values of a 0 might not be 

the true permanent Indentations. They were the values which 
could make the power law given by Equation (3-4) fit the 
total data under the unloading path. In fact, the load 
correspond I ng to the value of a cr = 3.16x10~ 3 In. obtained 
In [14] Is about 200 lb. for 0.5 In. Indenter, which Is 
apparently too high. The value of a p = 6.564x10" 4 In. 

obtained In this study, which corresponds to about 20 1 b of 
loading, seems more reasonable as a critical value In 

Indentation. For comparison, the relations between 

unloading rigidity s and maximum Indentation a m using 
Equation (3-7) with a. cr - 3.16x10" 3 In. and Equation (3-8) 

' with a p = 6.564x10 -4 In., respectively, are plotted In 

Figure 3,19. It Is Interesting to see that these two 
equations give almost the same values of s up to a m - 4x10” 3 
in. which Is approximately the maximum indentation before 
failure could occur to the specimen. The advantage of using 
Equation (3-7) for the formulation of the unloading law is 


f 
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Figure 3.19 Unloading rigidity s as function of maximum 
i ndentat i on 
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that the value of s Is constant for any a m once the the 
Indentation passes a cr> and only one unloading test Is 
necessary to determine a cr provided the load Is high enough 
to produce permanent Indentations. The use of Equation (3- 
9) needs performing many tests to obtain a proper relation 
between a 0 and a m according to Equation (3-8). However, It 
should be noted that Equation (3-7) Is valid only If q = 5/2 
Is used In the unloading equation (3-4), while Equation (3- 
8) has no such restriction. 
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CHAPTER 4 
IMPACT EXPERIMENTS 

High velocity Impacts usually result In very small 
contact time and the material under Impact loadings may 
behave differently from static contact due to the strain 
rate effect. The statically determined contact laws 
presented In the previous chapter thus must be verified 
experimentally before it can be appl led to the Impact 
analysis. Wang [15] has conducted many Impact experiments 
on laminated composite beams and plates using spherical 
steel balls as Impacters. The strain response histories at 
various points on the specimens were recorded and compared 
with the finite element analysis with which the contact laws 
obtained by Yang [14] was l ncorporated. The results showed 
that the test data agreed with the predictions using the 
statical Indentation laws quite well. In this chapter, an 
attempt was made to measure the contact force directly so 
that the applicability of statical contact laws In Impact 
analysis can be further evaluated. 
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4.1 Experimental Procedure 

A 6 In. by 4 In. laminated plate cut from a [0°/45°/0°/“ 
45 0 /O°] 2 s graphite/epoxy panel was used as the Impact 
target. The 0 0 -d! recti on was arranged to parallel the long 
side of the plate. Seven strain gages (Micro Measurement 
Company TYPE EA-1 3-062 AQ 350) were placed at different 
locations as shown In Figure 4.1 to record the dynamic 
strain histories. One of the gages was placed on the 
surface directly opposite to the Impact point to trigger the 
oscilloscope. This plate was hung with two strings at tv/o 
corners to achieve the free boundary condition. 

The projectile was made of an Impact-force transducer 
with a spherical steel cap of 0.75 Inch In diameter glued on 
the Impact side and a steel rod of 5/8 Inch In diameter 
glued on the other side as shown In Figure 4.2. It was then 
attached to a thin rod to form a pendulum which could 
produce Impact velocities up to 150 In/sec. The total mass 
of the projectile Is 0.000181 lb-sec 2 / In . 

The schematic diagram for this Impact experimental set-up 
is shown In Figure 4.3. Signals from gages and transducer 
were amplified by a 3A9 Textron I x amplifier and displayed on 
the screen of an oscilloscope. 
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Figure 4.1 Laminate dimension and strain gage locations 





0.75 in Did. 
Spherical Cap 


Transducer, 


|“ d,?'7 



5/8 in Dia. 
Steel Rod 




(a) Impact-Force Transducer 


(b) Projectile 


Figure 4.2 Graphical illustration of Impact projectile 
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Figure 4.3 Schematical diagram for the Impact 
experimental set-up 
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4.2 Calibration of Impact-Force Transducer P/IQK rJ 

OF POOR QUAUTV 

The Impact-force transducer used was Modal 200A05 
marketed by PCB Piezotronlcs Inc. Some of It's 
specifications are shown In Table 4.1 [30], The structure 
of this transducer contains two thin quartz disks operating 
In a thickness compression mode and sandwiched between 
hardened steel cylindrical members. A built- In amplifier 
can reduce the high Impedance of the voltage from the quartz 
element and provides an output voltage which can be read out 
on oscilloscope, recorder, etc,. The Impact force Is then 
computed using the equation, 

F “ Vp/cp (4-1 ) 

where V F Is the output voltage and c F Is the sensitivity of 
the transducer. Since the value of c F In Table 4.1 was 
obtained under quasi -static condition [30], It must be 
verified under Impact condition first so that later the 
results from Impact experiment can be correctly Interpreted. 

A circular cylindrical steel rod of 2 Inch In diameter 
and 1.19 Inch long hung on strings was used as the Impact 
target to calibrate the transducer. The acceleration of the 
rod was measured by using a Model 302A accel erometer which 
was mounted on the end of the rod opposite to the Impacted 
end as shown In Figure 4.4. The total weight of the target 
Is 1 ,105 lb. 
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Table 4. 1 

Specifications for Model 200A05 Impact-Force Transducer 


Range, Compression 
(5V output) 

lb. 

5,000 

Maximum Compression 

lb. 

10,000 

Resolution (200 /zV p-p noise) 

1b. 

0.2 

St 1 ffness 

Ib/Mln 

100 

Sens 1 1 1 v I ty 

mV/ lb 

1 .0 

Resonant Frequency 
(no load) 

Hz 

70 , 000 

Rise Time 

/i sec 

10 

Discharge Time Constant 
(T.C.) 

sec 

2,000 

Low-Frequency (-5%) 

Hz 

0.0003 

Linear I ty ,B.F.S,L. 

% 

1 

Output Impedance 

ohms 

100 

Excitation (thru C.C. diode) 

VDC/mA 

+18 to 24/2 to 20 

Temperature Coefficient 

%/°F 

0.03 

Temperature Range 

Op 

-100 to +250 

Shock (no load) 

g 

10,000 
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Using Equation (4-1) and 

a “ V a / c 0 (4-2) 

F = ma (4-3) 

we obtain 

Cp « (c Q /m) (V F /V a ) (4-4) 

where V Q and c Q are the output voltage and the sensitivity 
of the accel erometer , respectively, a Is acceleration of the 
target, and m is the mass of the target. 

When Impacting a metal projectile on a metal target with 
no pad on the impact surface, a high frequency ringing can 
be seen at the output of the transducer. In order to obtain 
smooth output curves, a soft pad was placed on the impact 
region of the target to eliminate the high frequency 
ringing. The cause of this ringing phenomenon will be 
discussed later. Typical output voltages of transducer and 
accel erometer read from the oscilloscope are shown in Figure 
4.5. Values of V F were plotted vs the corresponding values 
of V 0 taken from these two curves at several discrete points 
In time and then fitted into a straight line as shown In 
Figure 4.6. The slope of this line represents the ratio of 
V F /V a which is then substituted in Equation (4-4) to 
calculate the sensitivity c F . 



(mV) 
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Assuming the sensitivity of the accelerometer c„ Is 
correct, and using Equation (4-4) and the test data, the 
average value of c F calculated was 0.494 mV/lb.. A 
comparison with the value of 1.0 mV/ lb from Table 4,1 shows 
that the test result has more than 50% error. However, 
since the quartz elements are located at the center of the 
projectile while the impact force Is applied at the end, we 
were not certain that the force ht story picked up by the 
quartz elements did represent the real history of the Impact 
force. The following simple analysis was performed to 
examine this uncertainty. 

Consider a 1 In. long steel rod with free-free boundary 
conditions. For a Impulse loading given by 

F(t) = F 0 EXPC-(t-r) a /4b 2 ) ] (4-5) 

at one end, the force history at the midpoint of the rod, 
F m (t), was computed and plotted In Figure 4,7 together with 
the applied force history. It should be noted that the 
values of Fq - 1000 lb,, r = 200x10“ 6 sec. and b 40x10“ 6 
sec. were chosen in Equation (4-S) so that the applied 
force history Is similar to the experimental loading 
hlstroy. From Figure 4*7, It can be seen that F m (t) Is only 
about half of the applied force F(t). The average ratio of 
F m (t)/F(t) was obtained to be 0.498, which Is very close to 
the value of c F obtained previously. The accelerations at 
the two ends and the midpoint of tho rod were also 
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calculated and plotted In Figure 4.8* It shows that the 
magnitudes of acceleration at any position of the rod have 
virtually no difference, This Indicates that the 
accel eromoter did measure the real acceleration of the 
target while the Impact-force transducer only picked up the 
force history at the point of it's own position. In other 
words, the wave motion In the projectile can not be 
neglected, hence It must be treated as an elastic body, 

Repeating the previous analysis by changing the Impulse 
loading of Equation (4-B) to 

F(t) «» F 0 s I n(7rt/b) (4-6) 

and letting F 0 * 1000 lb. and b « 400x10” 5 sec,, we obtain 
the force history at the midpoint of the rod as shown In 
Figure 4.9. Comparing Figure 4.9 with Figure 4.8, It Is 
clear that the Initial slope of the Impulse forcing function 
would affect the amplitude of ringing. The steeper the 
Initial slope Is, the higher the amplitude of ringing will 
be. When Impacting the steel projectile on graph 1 te/epoxy 
surface, thts ringing phenomenon was also observed. 
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4.3.1 Plate Finite Element 

A 9-node Isoparametric plate finite element (see Figure 
4.10) developed by Yang [31] based upon the laminate thaory 
of Whitney and Fagano [18] was used to model the dynamic 
motion of the laminated plate. At each node there are five 
degrees of freedom. Among them, u° , v° and w are 
displacement components of mid-plane In the x-,y- and z- 
direction, respectively, and <*> x and <P V are rotations of the 
cross-sections perpend l cul ar to the x- and y-axls, 
raen^t-Uiaiu Prvr. o i- r t r. laminates. th® flexural 

deformation Is uncoupled from the In-plane extenslonal and 
shear deformat Ions , and hence, the degrees of freedom 
corresponding to u° and v° can be neglected In the 
transverse Impact problem, 

The Isoparametric plate finite element Is developed using 
the following shape functions; 

For corner nodes; 

S,w(1/4) (1+$o)(1+7?o) C^q+ 770-1 ) + (1 /4)(1H? a ) (1-r? 2 ) (4-7) 

For nodes at $ = 0 and 77 = ±1 : 

S , - ( 1 /2) ( 1 -$ 2 ) (7?g+?? s ) (4-8) 


For nodes at (? - ±1 and 7 ? = 0; 
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(b) DISTORTED ELEMENT 


Figure 4.10 9-node Isoparametric plate element 
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For the center node: 

S ( -(1/2) (1~$ a ) C 1 — T7 54 ) (4-10) 

In the above shape functions, $ and r? are normalized 
local coordinates, and 

$a * » Vo « W0\ (4-11 ) 

where and 77 , are the natural coordinates of node 1 
(Figure 4,10), 

Using the shape functions, the plate dl spl acements w» 
and <P y are approximated by 

w 

**■- 2 [S,J{q p h (4-12) 

1-1 

where { q p > , Is the nodal displacement vector at node 1 and 
3x3 

[S] , = S| [13 (4-13) 

The stiffness and mass matrices are obtained by numerical 
integration using Gauss quadrature. Following standard 
finite element procedures, the system stiffness matrix [K p ] 
and mass matrix [M p ] are assembled from the element 
matrices. The equations of motion are expressed in matrix 
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form as 



CM p Hq p } + [K p Hq p } - (P p > 


(4-14) 


where 


(P p > T ■“ {0, •••,F, •••„<» (4-15) 

Is the force vector In which F Is the contact force 
associated with the degree of freedom correspond! ng to the 
w-dlspl acement at the Impact point. The subscript p In 
Equations (4-12) through (4-15) denotes those are quantities 
corresponding to laminated plate. 


4.3.2 Modeling of Projectile 

In Section 4.2 we showed that In order to Interpret the 
experimental transducer response, It Is necessary to treat 
the projectile as an elastic body. A higher order rod 
finite element developed by Yang and Sun [32] was used to 
model the projectile. This element has two degrees of 
freedom at each node, namely the axial displacement u and 
It's first derivative 9u/9x. It has been shown that this 
higher order element Is far more superior than the elements 
with less degrees of freedom In the analysis of dynamic 
problems. The displacement function Is taken as 


u = a t + a 2 x + a 3 x 2 + a 4 x 8 


(4-16) 
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where a, are constant coefficients. Solving these 
coefficients In terms of the nodal degrees of freedom and 
substituting Into Equation (4-16), we obtain 

u - {N)T{q P >. (4-17) 

where 

{q r }* T - (<u) 1f <0u/dx) 1# <U) 2 , (8u/0x) a > (4-18) 

Is the vector of element nodal degrees of freedom, and 

{N) T « { f -j (x) , f a (x), f a (x) , f*(x>) (4-19) 

In which 

f, (x) » (1 - x/L) 2 ( 1 + 2 x/L) 
f 2 (x) ■ x(1 - x/L ) 2 
f a (x) » X 2 /L 2 (3 - 2x/L) 
f 4 (x) « x a /L(x/L - 1 ) 

are shape functions. The subscript r In Equation (4-17) 
denotes quantities correspond! ng to the rod. 

Using variational principle, the equations of motion for 
one element are obtained as 

Cm r ]{q r } 0 + [k r ]{q r } e - (p r } 0 (4-20) 

where {p r } 0 Is the vector of the generalized forces 
associated with the nodal degrees of freedom {q r } 0 , Cm r 3 Is 
the element mass matrix whose entries are given by 
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(m r ) , j = pA^f,fjdx I,j - 1,2, 3, 4 (4-21) 

and [k r ] Is the element stiffness matrix whose entries are 
given by * 

(k r ) , j = EAj^f , ' f j 'dx I, J « 1,2, 3, 4 (4-22) 

In Equations (4-21) and (4-22), p, E and A are mass density, 
Young's modulus and cross-sectional area of the projectile, 
respectively, and L Is the length of the element. The 
explicit forms of [k r ] and [m r 3 are given by 


and 


[k r 3 


[m r 3 



' 36 

3L 

-36 

3L 

EA 

3L 

4L 2 

-3L 

-L 2 

30L 

-36 

-3L 

36 

-3L 


3L 

-L 2 

-3L 

4L 



' 156 

22L 

54 

-13L 

pAL 

22L 

4L 2 

13L 

-3L 2 

420 

54 

13L 

156 

-22L 


-13L 

-3L 2 

-22L 

4L 2 


(4-23) 


(4-24) 


Following the usual manner, the system stiffness and mass 
matrices are assembled from the element stiffness and mass 
matrices, and the system equations of motion are expressed 
as 


CM r ] {q P > + [K r ] { q r } = {P r > 

where 

{P r } T = (F,0, • • • ,0} 

In which F Is the contact force 
of the projectile. 
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(4-26) 
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4.4. Results and Discussion 

The 6 In. by 4 In. graph I te/epoxy laminate was modeled 
by’ 140 (14 x 10 mesh) plate elements while the projectile 
was modeled by 20 rod elements (see Figure 4.11). The two 
sets of equations (4-14) and (4-25) along with the contact 
laws given by Equations (3-1), (3-3) and (3-11) were solved 
simultaneously. The finite difference method with At = 0.2 

'iSec. was used to Integrate the time variable. A coarser 
finite element mesh for plate was used and It was found that 
the present mesh yielded converged solutions. A 3- 
D! mens tonal analysis using 112 ax I symmetric finite elements 
to model the projectile was also performed, and the results 
showed the the response at the midpoint of the projectile to 
have no significant difference comparing with the one 
obtained by using rod elements. 

An Impact velocity of 115 In/sec was used In the 


experiment . 
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gages. The results obtained using the finite element 
methods and the contact laws are also shown In these 
figures. It Is evident that the finite element solutions 
agree with the experimental data very well* 

In Figure 4.18, the experimental transducer responses and 
the computed transducer responses using finite element are 
plotted against time as curve I and curve II, respectively. 
The computed contact force history Is also plotted as curve 
III, It can be seen that the magnitudes of curve I and 
curve II agree fairly well. The frequencies of ringing for 
these two curves, however, are quite different. For the 
finite element results, the time interval between two 
consecutive peaks of ringing Is approximately equal to the 
time that the longitudinal stress wave needed to travel the 
distance between two ends of the projectile. This Indicates 
that the ringing Is simply caused by the transient wave 
travelling back and forth In the projectile. 

From Figure 4.18 we can see that curve I has exact 9 
peaks In 180 microseconds, and the time Interval between two 
consecutive peaks Is about 20 microseconds. It Is noted 
that this transducer has a rise time of 10 microseconds (see 
Table 4.1), which Is the time it needs to reach the maximum 
response. Any Input signal with period smaller than twice 
of this value will be smoothed out by the transducer, and 
the output signal may appear to have lower frequency. In 
other words, the period of the output signal will be at 
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Figure 4.12 Strain response history at gage No.1 
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Figure 4.13 Strain response history at gage No. 2 
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Figure 4.14 Strain response history at gage No. 3 
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least 20 microseconds. This might explain the lower 
frequency of ringing In the output voltage from the 
transducer . 

The total duration of contact for this Impact test Is 
about 800 microseconds, and multiple contact Is also 
observed from the test data. Figure 4.19 shows the 
experimental transducer responses and the computed 
transducer responses up to 800 microseconds. Although these 
two results do not matched very well after the end of the 
first contact, It Is evident that the finite element 
analysis does predict the multiple contact phenomenon, and 
the calculated total duration of contact Is also 
approximately the same as the test result. 

Figure 4.20 presents a number of deformed configurations 
of the laminated plate after Impact. It Is seen that at the 
point of Impact, there Is a strong discontinuity In slope of 
the transverse displacement Indicating the presence of a 
significant transverse shear deformation. 
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Figure 4.19 Transducer response histories from 

experimental and finite element results up 
to 800 microseconds 
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TIME s 10/XSEC. 


TIME = 60/xSEC, 
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Figure 4.20 Deformed configurations of laminated pi sate 
after impact 
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CHAPTER 5 

SUMMARY AND CONCLUSION 

The laminate theory developed by Whitney and Pagano was 
employed for studies of harmonic wave and propagation of 
wave front In a [O°/45 0 /O 0 /-45 0 /O 0 ] 2 s graphite/epoxy 
laminate. The dispersion properties of flexural waves were 
Investigated. The wave front surface was constructed using 
ray theory, It was shown that due to the anisotropic 
properties of composite laminate, the transient wave would 
propagate with different velocities In different directions. 
The growth and decay of the wave front strength were also 
discussed. 

The contact laws between 0.5 Inch and 0.75 Inch spherical 
steel Indenters and the graph! te/epoxy laminate were 
determined experimentally by means of a statical Indentation 
test. Loading, unloading and reloading curves were fitted 
Into power equations. Linear relation was found between the 
permanent Indentation and the maximum Indentation at 
unloading, which Is seen to be Independent of the size of 
Indenters, This relation was then used to determine the 
coefficient of the unloading law. It was demonstrated that 
there was no need to perform reloading experiments once the 
loading and unloading laws were established. Test results 
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showed loading and reloading curves followed the power laws 
with power Indices of 1.5 very well, while the power Indices 
for unloading curves varied from 1,5 to 2.5, 

The statically determined contact laws were Incorporated 
Into an existing 9-node Isoparametr lc plate finite element 
program to study the dynamic response of a graphite/epoxy 
laminated plate subjected to Impact of a hard object. An 
Impact experiment was conducted to verify the validity of 
statical contact laws In the dynamical Impact analysis. It 
was shown that the strain responses predicted using the 
finite element method agreed with the test results very 
well. The contact force history of the Impact test was 
measured by an Impact-force transducer, which was also seen 
to match the finite element result In magnitude as well as 
contact duration. 


The Indentation tests have been used ever since the 
beginning of the century to determine the static and dynamic 
hardnesses of metals In'terms of the applied loading, the 
size of the Indenter, and the chordal diameter of the 
permanent Indentation [33]. If similar systematic 
Indentation tests are performed on the laminated composite 
materials, then the relations between contact coefficients 
and the sizes of the Indenters could be determined more 
rigorously, and the usefulness of the contact laws could be 
further extended. 


ORIGINAL page is 

OR POOR QUALITY 


115 


As th© verification of the contact laws has been limited 
to low velocity Impacts In this study, their accuracy under 
high velocity Impact conditions Is not clear. Besides the 
contact behavior which may be significantly different from 
the static one, the damage Induced by waves could be quite 
extensive which needs to be Included In the analysis. While 
the present study tried to establish experimental ly contact 
laws which can be used In the analysis of low velocity 
Impact, the damage of laminate due to Impact loading has not 
been discussed, It Is apparent that more work needs to be 
done so that the failure mechanism In laminated composites 
due to Impact can be better understood. Stress waves 
propagating In thickness direction, which may be responsible 
for the del ami nation of laminates, Is one of the Important 
subjects that should be Investigated, Strength and fatigue 
life degradations of laminates after Impact, which have been 
examined briefly by Wang [15], also need more extensive 
study. 
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APPENDIX 

COMPUTER PROGRAM AND USER INSTRUCTIONS 

The computer program used In this research was wrttten 
following the program by Professor R, L. Taylor [34] with 
some necessary modification In order to solve the Impact 
problems of laminated plates, A brief Instruction of the 
Input data for solving the Impact problem specified In 
Chapter 4 of this report Is given In this apppendlx. The 
detailed descriptions of data Input as well as the macro 
instructions for solving various types of problems can be 
found In [34], The listing of Input Is shown at the end of 
this appendix, followed by the listing of program. 

I. Title and control Information: 

1. Title card-Format (20A4) 

Col umns Poser 1 pt 1 on 

1-4 Must contain FECM 

5- 80 Alphanumeric Information to be printed with 

output as page header, 

2. Control information card-Format (615) 

Col umns Descr I pt i on 

1-5 Number of nodes (NUMNP) 

6- 10 Number of elements (NUMEL) 


r 
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11-15 

Number of 

layers (LAYER) 0F P00R Q UAL,TY 

16-20 

Spatial dimension (NDM) 


21-25 

Number of 

unknowns per node 

(NDF) 

26-30 

Number of 

nodes per element 

(NEN) 


II. Mesh and Initial Informat ton: 

The Input of each segment In this part of data is 
controlled by the alphanumeric value of macros, which must 
be followed Immediately by the appropriate data. Except for 
the END card which must be the last card of this part, the 
data segemnts can be In any order. Each segment Is 

terminated with blank card(s). The meaning of each macro Is 
given by the following: 

Macro Data to be Input 

C00R Coordinate data 

ELEM Element data 

BOUN Boundary condition data 

MATE Material data 

ROD Initial condition of the projectile 

EXPE Experimental Indentation laws data 

END Must be tho last card of this part, terminates 

mesh and Initial Information input. 

1, Coordinate data-Format(2I5,2F10.0) 

Columns Description 

1-5 Nodal number 

6-10 Generation Increment 
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21-30 Y-coordlnate 

2. Element data-Format ( 1 1 1S) 

Co 1 umns Description 

1-5 Element number 
6-10 Node 1 number 
11-15 Node 2 number 
etc. 

46-50 Node 9 number 
51-55 Generation Increment 

3. Boundary condition data-Format (71'5) 

Columns Description 

1 -5 Node number 
6-10 Generation Increment 
11-15 DOF 1 boundary code 

1 6-20 DOF 2 boundary code 

21-25 DOF 3 boundary code 

26-30 DOF 4 boundary code 

31-35 DOF 5 boundary code 


4. Initial condition of the project 1 le-Format(2I5,F10.0) 
Col umns Descr I pt 1 on 

1-5 The node at which the projectile hits 


6-10 DOF correspond! ng to the direction of Impact 
11-20 Initial Impact velocity 
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5. Experimental Indentation laws data-Farmat(4F10,0) 
Columns Descriptio n 

1-10 Contact coefficient k 
11-20 Critical Indentation a p 

21-30 Constant s p of Equation 3-9 

31-40 Power Index q of the unloading lav/ 

6. Material data 

Card 1 -forma t(315,F10,0) 

Columns Descr Ipt Ion 

1-5 Order of Gauss quadrature for the numerical 
integration of the bending energy 
6-10 Order of Gauss quadrature for the numerical 
integration of the transverse shear energy 
11-15 Order of Gauss quadrature for strain outputs 
at Gauss points If >0 

at nodal points if <0 

16-25 Total thickness of the laminate 

Card 2~Format(7F10.0) 

Columns Descr I pt I on 
1-10 Mass density 
11-20 Poisson's ratio i/ 12 
21-30 Longitudinal Young's modulus E t 
31-40 Transverse Young's modulus E 2 
41-50 Shear modulus G 12 

11-20 Shear modulus G 13 

Shear modulus Ga S 


11-20 
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Card 3, 4, * • * Format* I5,F5.0,F10,0) 

Colurins Dese rt, j ot l on ORIGINAL PAQg KJ 

1-5 Layer number 0F QUALITY 

6-10 Fiber angle 
11-20 Thickness of the layer 

III. Macro instructions: 

The first Instruction must be a card with MACR In columns 
1 to 4. The macro Instructions needed to solve the problem 
specified In Chepter 4 of this report are shown In the 
listing of Input. Cards must be Input In the precise order. 
The following Is the explanation of each macro: 

Columns Columns Columns 


1-4 

5-10 

11-15 

Description 


LMAS 



Lumped mass formulation 


DT 


V 

Set time Increment to va 

lue V 

LOOP 


N 

Execute N times the Instructions 




between this macro and macro NEXT 

TIME 



Advance time by DT value 


RODP 


N 

Integration of the equat 

ions of 




motion using the finite 

d I f ference 




method. Contact force, I 

ndentat Ion 




and element strain will 

be stored 




stored every N steps in 

loop 

DISP 


N 

Nodal displacements will 

be stored 




every N steps In loop 


NEXT 



End of loop instructions 
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END End of macro program Instruct ions 

IV, Termination of program execution 

A card with STOP In columns 1 to 4 must be supplied nt 
the end of the Input data In order to properly terminate tho 
execution. 

The values of contact force, Indentation, element strain, 
nodal displacement and the response of the projectile at 
each requested output time step are stored In program files 
which can be saved (say, copy to a magnetic tape) at the end 
of execution. Three program files, !,©.* tapes, tap©8 and 
tapes are used for data saving? 

TapeS: Nodal displacement - Format (6E1 2,4) 

Nodal displacements, from node 1 to node NUMNP, are saved 
on tape3 at each requested output time step according to the 
format . 

TapeS: Element strain - Format(2I6,5E12,4) 

Element strains, from element 1 to element NUMEL, and 
then from node 1 to node NEN of each element, are saved on 
on tape© at each requested output time step. 

Columns Data saved 

. « ! ' '■ ■■■» m ' ■■« ! ■■■- 1 ■ |«> II « .1^ n 

1-6 Element number 

7-12 Node number of dement 

13-24 Bending strain 

Bending strain /c v 


25-36 
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37-48 Bending strain fc xy 

49-60 Transverse shearing strain i y * 

49-60 Transverse shearing strain *V Xis 

Tape9: Contact force, indentation and the response of the 
projectile - Format (6E1 2. 4) 

The following Information Is saved on tape9 at each 
requested output time step: 

Columns Data saved 
1-12 Contact force 
13-24 Indentation 

25-36 'Transducer' response (see Chapter 4) 

37-48 Displacement of the projectile at the Impacted end 

37-48 Velocity of the projectile at the Impacted end 

37-48 Acceleration of the projectile at the Impacted end 
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LISTING OF INPUT DATA 


FECM *#L0H UELOCITY IMPACT OF LAMINATED PLATE** 
G 09 HO 20 2 5 9 

COOR 


1 

1 

0.0 

0.0000 

7 

1 

1.5 

0.0000 

23 

1 

4.5 

0.0000 

29 

0 

6.0 

0.0000 

30 

1 

0.0 

0.2500 

36 

1 

1.5 

0,2500 

52 

1 

4.5 

0.2500 

50 

0 

6.0 

0.2500 

59 

1 

0,0 

0.5000 

65 

1 

1.5 

0.5000 

61 

1 

4.5 

0.5000 

87 

0 

G.O 

0.5000 

03 

1 

0.0 

0.6875 

94 

1 

1.5 

0.6875 

110 

1 

4.5 

0.6875 

116 

0 

G.O 

0,6875 

117 

J 

0.0 

0.0750 

123 

I 

1.5 

0.8750 

139 

1 

4.5 

0,0750 

145 

0 

6.0 

0.8750 

146 

1 

0.0 

1.0625 

152 

1 

1.5 

i ntsac 

x * vuuu 

163 

1 

4.5 

1.0625 

174 

0 

6.0 

1.0625 

175 

1 

0.0 

1.2500 

181 

1 

1.5 

1,2500 

197 

1 

4.5 

1.2500 

203 

0 

G.O 

1.2500 

204 

1 

0.0 

1.4375 

210 

1 

1.5 

1.4375 

226 

1 

4.5 

1.4375 

232 

0 

6.0 

1.4375 

233 

1 

0.0 

1.6250 

239 

1 

1.5 

1.6250 

255 

1 

4.5 

1.6250 

261 

0 

6.0 

1.6250 

262 

1 

0.0 

1.8125 

268 

1 

1.5 

1,8125 

284 

1 

4.5 

1.8125 

290 

0 

G.O 

1.8125 

291 

1 

0.0 

2.0000 

297 

1 

1.5 

2.0000 

313 

1 

4.5 

2.0000 

319 

0 

6.0 

2.0000 

320 

1 

0.0 

2.1875 

326 

1 

1.5 

2.1875 

342 

1 

4.5 

2.1875 

348 

0 

6.0 

2.1875 

349 

1 

0.0 

2.3750 

355 

1 

1.5 

2.3750 

371 

1 

4.5 

2.3750 

377 

0 

G.O 

2.3750 

378 

1 

0.0 

2.5625 

384 

1 

1.5 

2.5625 

400 

1 

4.5 

2.5625 

406 

0 

0.0 

2.5625 

407 

1 

0.0 

2.7500 

413 

1 

1.5 

2.7500 

429 

1 

4.5 

2.7500 

435 

0 

6.0 

2.7500 

436 

1 

0.0 

2.9375 

442 

1 

1.5 

2.9375 
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458 

1 

4.5 

2.9375 

4G4 

0 

6.0 

2.9375 

465 

1 

0.0 

3.1250 

471 

1 

1.5 

3.1250 

487 

1 

4.5 

3.1250 

493 

0 

6.0 

3.1250 

494 

1 

0.0 

3.3125 

500 

1 

1.5 

3.3125 

51G 

1 

4.5 

3.3125 

522 

0 

6.0 

3.3125 

523 

1 

0.0 

3.5000 

529 

1 

1.5 

3.5000 

545 

1 

4.5 

3.5000 

551 

0 

6.0 

3.5000 

552 

1 

0.0 

3.7500 

558 

1 

1.5 

3.7500 

574 

1 

4.5 

3.7500 

580 

0 

6.0 

3.7500 

581 

1 

0.0 

4.0000 

587 

1 

1.5 

4.0000 

603 

1 

4.5 

4.0000 

609 

0 

6.0 

4.0000 

ELEN 

1 

1 

3 61 

59 2 

15 

59 

61 119 

117 60 

29 

117 

119 177 

175 118 

43 

175 

177 235 

233 176 

57 

233 

235 293 

291 234 

71 

291 

293 351 

349 292 

85 

349 

351 409 

407 350 

99 

407 

409 467 

465 408 

113 

465 

467 525 

523 466 

127 

523 

525 583 

581 524 

BOUN 

1 

1 

-I -1 

0 0 

609 

0 

1 1 

0 0 

ROD 

305 

3 

115.0 


EXPE 

1912000. 

0.0006564 

0.094 

NATE 

3 

3 

-3 

.106 

0.000148 

0.3 

17500000. 

1 

0 . 

0.0053 


2 

45. 

0.0053 


3 

0 . 

0.0053 


4 

-45. 

0.0053 


5 

0 . 

0.0053 


6 

0 . 

0.0053 


7 

45. 

0.0053 


8 

0 . 

0.0053 


9 

-45. 

0.0053 


10 

0 . 

0.0053 


11 

0 . 

0.0053 


12 

-45. 

0.0053 


13 

0 . 

0.0053 


14 

45. 

0.0053 


15 

0 . 

0.0053 


16 

0 . 

0.0053 


17 

-45. 

0.0053 


18 

0 . 

0.0053 


13 

45. 

0.0053 


20 

0 . 

0.0053 



32 

60 

30 

31 

2 

90 

118 

88 

89 

2 

148 

176 

146 

147 

2 

206 

234 

204 

205 

2 

264 

292 

262 

263 

2 

322 

350 

320 

321 

2 

380 

408 

378 

379 

2 

438 

466 

436 

437 

2 

496 

524 

494 

495 

2 

554 

582 

552 

553 

2 


o 

o 


a.o 


1150000. 800000. 800000. 


800000. 


END 
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MACR 

LMAS 

DT 

LOOP 

TIME 

RODP 

DISP 

NEXT 

END 

STOP 


.SE-S 

10 


5 

5 


ORIGINAL FA&2 
OF POOR QUALITY 



ORIGINAL PAGE 
OF POOR QUALITY 


129 


LISTING OF PROGRAM 


PROGRAM MAINdNPUT, OUTPUT* TAPE5*INPUT, TAPE6»0UTPUT* TAPE2* TAPE3* 
1 TAPE8, TAPES) 

Off#*-* MAIM PROGRAM 
LOGICAL PCOMP 
COMMON /PRSIZE/ MAX 

COMMON /CTDATA/ 0 , HE AD ( 2 0 ) , NUMNP » NUMEL * LA YER , NEQ » IPR 
COMMON /'LABELS/ PDIS(6),A(G),BC(2),DI(G),CD(3),FD(3) 

COMMON /LODATA/ NDF * NDM, NEN* NST , NKM 
COMMON /PARATS/ NPAR(14),NEND 
DIMENSION TITL(20),WD(3) 

COMMON G(39000) 

DIMENSION M(39000) 

EQUIUALENCE (G(l).tl(l)) 

MAX=39000 
WD( 1 )=4HFECM 
WD(2)=4HMACR 
WD(3)=4HST0P 
999 READ <5* 1000) TITL 

IF(PCQMP(TITL(1)*WD(1) )) GO TO 100 
IF(PC0MP(TITL(1)*WD(2))) GO TO 200 
IF(PCOMP(TITL( 1) , WD(3) ) ) STOP 
GO TO 999 

100 DO 101 1=1*20 

101 HEAD(I)=TITL(I) 

R£AD(5» 1001 ) NUMNP* NUMEL. LAYER, NDM, NDF, NEN 
WRITE( G, 2000 ) HEAD, NUMNP, NUMEL, LAYER, NDM, NDF, NEN 
PDIS(2)=A(NDM) 

NST=NEN*NDF 
DO 110 1=1,14 
110 NPAR( I)=l 
NPAR(1)=1 

NPAR ( 2 ) =NPAR ( 1 ) +3*NST*IPR 
NPAR ( 3 ) =NPAR ( 2 ) +NDM*NEN*IPR 
NPAR ( 4 ) =NPAR ( 3 ) +NST 
NPAR ( 5 ) =NPAR ( 4 ) +NST« I PR 
NPAR ( G ) =NPAR ( 5 ) +NEN«NUMEL 
NPAR ( 7 ) =NPAR ( G ) +NDF*NUMNP 
NPAR ( 8 ) =NPAR ( 7 ) +NDM«NUMNP*IPR 
NPAR ( 9 ) =NPAR C 8 ) +NDF*NUMNP*IPR 
NPAR (10) =NPAR ( 9 ) +NDF*NUMNP 
CALL SETMEM(NPAR(9) ) 

CALL PZERO ( G ( 1 ) , NPAR ( 9 ) ) 

CALL PMESH ( M ( NPAR ( 3 ) ) , G ( NPAR ( 2 ) ) ,M(NPAR(5) ),MCNPAR(G) ), 

1 G ( NPAR ( 7 ) ) , G ( NPAR ( 8 ) ) , M ( NPAR ( 9 ) ) , NDF , NDM , NEN , NKM ) 

NPAR ( 1 0 )=NPAR ( 9 ) +NEQ 

NPAR (11) =NPAR (10) +NDF*NUMNP*IPR 

MEND=NPAR (11) +NEQ* IPR 

NE=NEND 

CALL SETMEH(NE) 

CALL PZERO ( G ( NPAR (10)), NE-NPAR (10)) 

GO TO 999 

200 CALL PMACR(G(NPAR(1)),G(NPAR(2)),M(NPAR(3)),G(NPAR(4)), 

1 M(NPAR(5) ) ,M(NPAR(G) ) »G(NPAR(7) ) »G(NPAR(8) ) ,M(NPAR(9) )» 

2 G(NPAR( 10) ) ,G(NPAR ( 11 ) ) , G(NE) , NDF, NDM, NEN, NST) 

CALL PZERO (G, MAX) 

GO TO 999 

1000 FORMAT ( 20 A4) 

1001 F0RMATUGI5) 

2000 FORMAT ( lHl * 20A4// 

1 5X,rC ONTROL INFORMATION S*// 


2 10X, 35HNUMBER OF NODAL POINTS =,16/ 

3 10X.35HNUMBER OF ELEMENTS =,IG/ 

4 10X.35HNUMBER OF MATERIAL LAYERS =, 16/ 

5 10X, 35HDIMENSI0N OF COORDINATE SPACE =, IG/ 


G 10X* 35HDEGREES OF FREEDOM FOR EACH NODE =,IG/ 
7 10X, 35HN0DES PER ELEMENT (MAXIMUM) =, IC) 

END 


MAIN 1 
MAIN 2 
MAIN 3 
MAIN 4 
MAIN 5 
MAIN G 
MAIN 7 
MAIN 8 
MAIN 9 
MAIN 10 
MAIN 11 
MAIN 12 
MAIN 13 
MAIN 14 
MAIN 15 
MAIN 1G 
MAIN 17 
MAIN 18 
MAIN 19 
MAIN 20 
MAIN 21 
MAIN 22 
MAIN 23 
MAIN 24 
MAIN 25 
MAIN 2G 
MAIN 27 
MAIN 28 
MAIN 29 
MAIN 30 
MAIN 31 
MAIN 32 
MAIN 33 
MAIN 34 
MAIN 35 
MAIN 3G 
MAIN 37 
MAIN 38 
MAIN 39 
MAIN 40 
MAIN 41 
MAIN 42 
MAIN 43 
MAIN 44 
MAIN 45 
MAIN 46 
MAIN 47 
MAIN 48 
MAIN 43 
MAIN 50 
MAIN 51 
MAIN 52 
MAIN 53 
MAIN 54 
MAIN 55 
MAIN 56 
MAIN 57 
MAIN 58 
MAIN 59 
MAIN SO 
MAIN 61 
MAIN G2 
MAIN G3 
MAIN G4 
MAIN G5 
MAIN GG 
MAIN G7 
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BLOCK DATA 
C***# BLOCK DATA 

COMMON /CTDATA/ 0* HEAD C SO >,NUMNP»NUMEL, LAYER* NEQ, IPR 
COMMDN /LABELS/ PDISC6), ACG)*BC(2),DI(6),CDC3'),FD(3) 

DATA 0/1H1/. IPR/1/ 

DATA PDIS/4H(I10,2H» ,4HF13.,4H4, ,4HGE13,4H.4) / 

DATA A/2H* 1 , 2H, 2, 2H, 3, 2H, 4, 2H, 5, 2H» S/ 

DATA BC/4H B.CiSH. / 

DATA DI/4H DIS,2HPL,4H UEL,2H0C,4H ACC,SHEL/ 

DATA CD/4H COO, 4HRDIN, 4HATES/ 

DATA FD/4H FOR, 4HCE/B, 4HISPL/ 

END 

C 

SUBROUTINE PMACR(IJL,XL,LD,P, IX, ID»X»F» JDIAG,DR,B,CT,NDF»NDM, 
* NEN»NST) 

MACRO INSTRUCTION ROUTINE 
LOGICAL Pu'MP 
COMMON G( 1 ) 

DIMENSION MCI) 

EQUIUALENCE (GC1),M(D) 

COMMON /CTDATA/ 0, HEAD(20),NUMNP,NUMEL> LAYER, NEQ, IPR 

COMMON /PROLOD/ PROP 

COMMON /THDATA/ TIME, DT, DDT, FORCE, ALPHA 

COMMON /ISWIDX/ TSU 

COMMON /PARATS/ NPARC 14) , NEND 

COMMON /RODATA/ UR,IQ,NDS 

DIMENSION UL(l),XL<l),LDCl),P(l),lXaMD(l),X(l),FU)# 

A JBIAfUl), DR(1)?B(1) 

DIMENSION WD(9)»CT(4» 1G)»LUE(9) 

DATA f 4D/4HL00P, 4HNEXT , 4HDT ,4HPR0P,4HLMAS,4HR0DP, 

1 4HSTRE, 4HDISP, 4HCHEC/ 

DATA NWD/9/,ENDM/4HEND / 

C.... INITIALIZATION 
DT = 0.0 
PROP = 1.0 
TIME = 0.0 
NNEQ = NDF«NUMNP 
NPLD = 0 
FORCE= 0. 

ALPHA- 0 1 

URITE(G,2001) 0,HEAD 
LL = 1 
LMAX = 1G 

CALL SETMEM ( NEND+LMAX*4* IPR ) 

CT(l.l) = WDC1) 

CTO, 1 ) = 1.0 
100 LL = LL + 1 

IFCLL.LT. LMAX) GO TO 110 

LMAX = LMAX + 1G 

CALL SETMEM <NEND+LMAX*4« IPR) 

110 READC5, 1000) (CT(J,LL), J“l,4) 

WRITECG, 2000) CCT(J,LL), J=l,4) 

IF C . NOT . PCOMP ( CT ( 1 , LL ) , ENDM ) ) GO TO 100 
CT(1,LL) = WDC2) 

NEND = MEND +LMAX#4«IPR 
LX = LL - 1 
DO 230 L=l, LX 

IFC.N0T.PC0MP(CTC1,L),WDC1))) GO TO 230 
J = 1 
K = L + 1 
DO 210 I=K, LL 

IFCPCOMPCCTC 1, I ) ,WD(1))) J = J + 1 
IFCJ .GT. 9) GO TO 401 
IFCPCOMPCCTC 1, I ) » WDC2) ) ) J = J - 1 
210 IFCJ.EQcO) GO TO 220 
GO TO 400 
220 CT (4,1) = L 
CT(4,L) = I 
230 CONTINUE 
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BLOC 1 
BLOC 8 
BLOC 3 
BLOC 4 
BLOC 5 
BLOC G 
BLOC r 
BLOC fJ 
BLOC 3 
I LOG 10 
BLOC 11 
BLOC 18 

nine i 
PMAC 2 
PMAC 3 

nine 4 

PMAC U 
PMAC G 
PMAC 7 
PMAC 0 
PMAC 0 
PMAC 10 
PMAC 11 
PMAC 18 
PMAC 13 
PMAC 14 
PMAC 1U 
PMAC 1G 
PMAC 1? 
PMAC 10 
PMAC 19 
PMAC 20 
PMAC 21 
PMAC 22 
PMAC 23 
PMAC 24 
PMAC 25 
PMAC 2G 
PMAC 2? 
PMAC 28 
PMAC 29 
PMAC 30 
PMAC 31 
PMAC 32 
PMAC 33 
PMAC 34 
PMAC 35 
PMAC 3G 
PMAC 3? 
PMAC 38 
PMAC 39 
PMAC 40 
PMAC 41 
PMAC 42 
PMAC 43 
PMAC 44 
PMAC 45 
PMAC 4G 
PMAC 4 7 
PMAC 48 
PMAC 49 
PMAC r,o 
PMAC 51 
PI'IAG ML- 
PHAC 53 
PMAC 51 
PMAC 55 
PMAC 5G 
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J» 0 

DO 840 L*1*LL 

IF(PC0MP(CTU»L),WD(1))) J « J + 1 
540 IF<PCONP(CT(l,L>,WDC5))) J * J - 1 
IFCJ„NE.O) GO TO 400 
LU a 0 
Ln 1 

599 DO 300 JaliNUD 

300 IPCPGOMP(CTQ,L),HD(J))) GO TO 310 
GO TO 330 
310 I n L - 1 

GO TO Ci,5,3,4»5»G,7,0,9),J 
C.... SET LOOP START INDICATORS 
1 LU « LU + 1 
LX b CT(4»L) 

LUECLU) n LX 
CTO, LX) a I. 

GO TO 330 

C.... LOOP TERMINATOR CONTROL 
5 N » CTC4.L) 

CTC3, L) « CTC3,L) + 1.0 
IF(CT(3,L),GT,CTC3»N)) LU = LU - 1 
IF(CT(3»L),LE.CT(3»N)) L * N 
GO TO 330 

C.,.. SET TINE INCREMENT 
3DTn CT(3»L) 

DDTb DT«DT 
GO TO 330 

C,... INPUT PROPORTIONAL LOAD TABLE 

4 NPLD « CTC3.L) 

PROP n PROPLDCO. » NPLD) 

GO TO 330 

C..,. FORM LUMPED MASS MATRIX 

5 ISM*»3 
CALL KMLIB 
GO TO 330 

C IMPACT 

G MDS=CT(3,L) 

IFCNDS.EQ.O) MDS«1 
CALL RODIPCT 
GO TO 330 

C.... PRINT STRESS/STRAIN UALUE 
7 ISM"4 
LX b LUECLU) 

IF ( At tOD ( QT ( 3 » LX ) , ANAX 1 < CT ( 3 , L ) , 1 . ) ) ) 330 , 71 , 330 
71 CALL FSTRLACUL, XL, LD, P, IX, ID, X, F, JDIAG, DR, B, NDF, NDM, NEN, NST, NNEQ) 
GO TO 330 

C. ... PRINT DISPLACEMENTS 
0 LX n LUE(LU) 

ir(AtlOD(CT(3,LX) , AMAX1(CT(3,L), 1, ) ) ) 330,81,330 
01 CALL FRTDISCUL, ID, Xv D, F, DR, NDM, NDF) 

GO TO 330 
C.,.. CHECK 

9 NRITECG, GOOD NEND, JDIAGCNEQ) 

RETURN 
330 u=u+i 

IKL.GT.LL) RETURN 
Gt) TO 893 

C.... PRINT ERROR FORMATS 

400 RRITE(G»4QQ0) 

RETURN 

401 MRITECB, 4001) 

RETURN 

C.... INPUT^OUTPUT FORMATS 
1000 FORMAT CfY1» IX, A4, 1X»2F5.Q) 

2000 FORMAT! 1 OX, A4, IX, A4, 1X,2G15.5) 

5001 FORMAT ( A1 , 20A4/V, 5X, 1GHMACRO INSTRUCTIONS/^, 15HMACR0 STATEMENT 
a , SX, lOHUARl’ABLE 1,5X, 10HUARIABLE 2) 

4000 FORMAT CSX, 4GH««PMACR ERROR 01»» UNBALANCED LOOP. NEXT MACROS ) 

4001 FORMAT C5X,45H»«PHACR ERRQE G2«» LOOPS NESTED DEEPER THAN 8) 


PMAC 57 
PMAC 58 
PMAC 53 
PMAC GO 
PMAC SI 
PMAC G2 
PMAC S3 
PMAC B4 
PMAC G5 
PMAC SR 
PMAC 6? 
PMAC GO 
PMAC S3 
PMAC 70 
PMAC 71 
PMAC 75 
PMAC 73 
PMAC 74 
PMAC 75 
PMAC 7G 
PMAC 77 
PMAC 78 
PMAC 79 
PMAC 80 
PMAC 81 
PMAC 85 
PMAC 83 
PMAC 84 
PMAC 85 
PMAC 8G 
PMAC 87 
PMAC 88 
PMAC 89 
PMAC 90 
PMAC 91 
PMAC 95 
PMAC 93 
PMAC 94 
PMAC 95 
PMAC 9G 
PMAC 97 
PMAC 98 
PMAC 99 
PMAC100 
PMAC101 
PMAC1QH 
PMAC103 
PMAC 104 
PMAC105 
PMAC10G 
PMAC 107 
PMAC 108 
PMAC109 
PMAC 110 
PM AC 111 
PMAC112 
PM AC 113 
PMAC 114 
PMAC115 
PMAC11G 
PNAC117 
PMAC! 18 
PMAC119 
PMAC1S0 
PMAC 121 
PMAC12S 
PMAC123 
PMAC 124 
PMAC1S5 
PNAC1HB 
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5001 FORMAT < 1H1 * ////5X* 32HCHECK MESH DATA AMD MEMORY SPACE// 
a 10X.12H MEMD «,I10//10X,12HJDIAG(NEQ) «,I10> 

EMD 

C 

SUBROUTINE P2ERO<U*NN) 

C***« ZERO REAL ARRAY 
DIMENSION UCNN) 

DO 100 N-liNN 
100 UCN) * 0.0 
RETURN 
END 
C 

SUBROUTINE SETMEMCJ) 

Caw# MONITOR AUAIABLE MEMORY IN BLANK COMMON 
COMMON /PRSIZE/ MAX 
K 6 J 

IF(K.LE.MAX) RETURN 
WRITECB* 1000) Kf MAX 
STOP 

1000 FORMAT CSX* 4SH##SETMEM ERROR 01** INSUFFICIENT STORAGE ,IN BLANK* 
* OH COMMON //17X.11HREQUIRED «, I0/17X* 11HAUAILABLE *, 18) 

END 

C 

LOGICAL FUNCTION PCOMPCA*B) 

C*#™ LOGICAL COMPARISON 
IF(A-B) 10,20*10 
10 PCOMP = .FALSE. 

RETURN 

20 PCOMP «• .TRUE. 

RETURN 

END 

C 

SUBROUTINE ACTCOLCA, B, JDIAG, NEQ, AFAC, BACK* ISS) 

C*n#„ ACTIUE COLUMN PROFILE SYMMETRIC EQUATION SOLUER 
LOGICAL AFAC* BACK, FLAG 
DIMENSION A(i)*B(l), JDIAG(l) 

C.... FACTOR A TO UT«D«U, REDUCE B 
FLAGa. FALSE. 

JR » 0 

DO 600 Jal.NEQ 
jn = JDIAGCJ) 

JH = JD - JR 
IS « J - JH + a 
IFCJH-2) 600*300* 100 
100 IFC. NOT. AFAC) GO TO 500 
IE « J - 1 
K = JR + 2 
ID « JDIAGCIS-1) 

C..., REDUCE ALL EQUATIONS EXCEPT DIAGONAL 

DO 200 I»IS* IE 
IR » ID 
ID » JDIAGCI) 

IH » NINO ( ID— IR— 1 » I-IS+1 ) 

IFCIH.GT.O) ACK)=ACK)-DOTCACK-IH)* ACID-IH), IH) 

200 K « K + 1 

C.... REDUCE DIGONAL TERM 

300 IFC. NOT. AFAC) GO TO 500 
IR = JR + 1 

IE = JD > 1 
K a j - JD 
DO 400 I a IR* IE 
ID a JDIAGCK+I) 

IFCACID) ) 301,400,301 

301 D a A(i) 

ACI ) = ACD/ACID) 

AC JD) a PKJD) - D*ACI) 

400 CONTINUE 

IFCACJD) )450*450,500 
450 IFCISS.NE. 0) GO TO 500 
IFCFLAG) GO TO 4S5 


PMAC1G7 

PMAC128 

PMAC1S3 

PEER 1 
PEER 2 
PEER 3 
PEER 4 
PEER S 
PEER G 
PEER 2 

SETM 1 
SETM a 
SETM 3 
SETM 4 
SETM 5 
SETM 0 
SETM ? 
SETM 0 
SETM H 
SETM 10 

PCOM 1 
PCOM 2 
PCOM 3 
PCOM 4 
PCOM 5 
PCOM 0 

peon ? 

PCOM U 

ACT!- 1 
ACTC 2 
ACTC 3 
ACTC 4 
ACTC S 
ACTC R 
ACTC ? 
ACTC 8 
ACTC 9 
ACTC 10 
ACTC 11 
ACTC 12 
ACTC 13 
ACTC 14 
ACTC 15 
ACTC 16 
ACTC 1? 
ACTC 18 
ACTC 19 
ACTC 20 
ACTC 21 
ACTC 22 
ACTC 23 
ACTC 24 
ACTC 25 
ACTC 26 
ACTC 27 
ACTC 28 
ACTC 29 
ACTC 30 
ACTC 31 
ACTC 32 
ACTC 33 
ACTC 31 

ACTC an 

ACTC 38 
ACTC 37 
ACTC 38 
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WRITEC6,4G0> 

ACTC 

39 

4G0 

FORMAT (//50H*#ACTCOL ERROR 01#* STIFFNESS MATRIX MOT POSITIUE » 

ACTC 

40 

1 8HDEFINITE) 

ACTC 

41 


FLAG=.TRUE. 

ACTC 

42 

•1GS 

WRITE(G,466) J,A(JD) 

ACTC 

43 

406 FORMAT < 32H NQNPOSITIUE PIUOT FOR EQUATION , I4.5X.7HP0UIT a. 

ACTC 

44 

- EGO. 10) 

ACTC 

45 

C * . . • 

REDUCE RHS 

ACTC 

4G 

500 

IF(BACK) B(J) = B(J) - DOTCA( JR+1 ) > B ( IS~1 ) » JH“1 ) 

ACTC 

47 

GOO 

JR = JD 

ACTC 

48 


IF(FLAG) STOP 

ACTC 

49 


IF (.MOT. BACK) RETURN 

ACTC 

50 

C< • • « 

DIUIDED BY DIAGONAL PIUOTS 

ACTC 

51 


DO 700 Isl.MEQ 

ACTC 

52 


ID = JDIAGCI) 

ACTC 

53 


IF (ACID) ) G50»700f G50 

ACTC 

54 

G50 

B( I ) = B(I)/A(ID) 

ACTC 

55 

700 

CONTINUE 

ACTC 

56 

C • . . » 

BACK SUBSTITUTE 

ACTC 

57 


J a NEO 

ACTC 

58 


JD = JDIAG(J) 

ACTC 

59 

800 

D = B(J) 

ACTC 

GO 


J = J - 1 

ACTC 

61 


IF(J.LE.O) RETURN 

ACTC 

62 


JR = JDIAGCJ) 

ACTC 

63 


IF( JD-JR.LE. 1 ) GO TO 1000 

ACTC 

G4 


IS = J - JD + JR + a 

ACTC 

G5 


K = JR - IS + 1 

ACTC 

GG 


DO 300 I=IS,J 

ACTC 

G7 

900 

B(I) = B(I) - A(I+K)*D 

ACTC 

68 

1000 

JD = JR 

ACTC 

69 


GO TO 800 

ACTC 

70 

p 

END 

ACTC 

71 


SUBROUTINE ADDSTF C A . S»P» JDIAG» LD»NST* NEL, FLG) 

ADDS 

1 

C##*# 

ASSEMBLE GLOBAL ARRAYS 

ADDS 

2 


LOGICAL FLG 

ADDS 

3 


DIMENSION A(1)»S(NST» 1)»P(1)» JDIAGCI )»LD(1) 

ADDS 

4 


DO 200 J=lf NEL 

ADDS 

5 


K = LD(J) 

ADDS 

6 


IFCK.EQ.O) GO TO 200 

ADDS 

7 


IF (FLG) GO TO 50 

ADDS 

8 


A(K)=A(K)+P(J) 

ADDS 

9 


GO TO 200 

ADDS 

10 

50 

L = JDIAG(K) - K 

ADDS 

11 


DO 100 1=1 » NEL 

ADDS 

12 


M = LD(I) 

ADDS 

13 


IFCM.GT.K .OR. M.EQ.O) GO TO 100 

ADDS 

14 


M = L + 11 

ADDS 

15 


A(N)=A(M)+S( I* J) 

ADDS 

16 

100 

CONTINUE 

ADDS 

17 

GOO 

CONTINUE 

ADDS 

18 


RETURN 

ADDS 

19 

p 

END 

ADDS 

20 

u 

FUNCTION DOTCA,B,N) 

DOT 

1 

Cwswtt 

UECTOR DOT PRODUCT 

DOT 

2 


DIMENSION A(l)»BCi) 

DOT 

3 


DOT = 0.0 

DOT 

4 


DO 100 I=l f N 

DOT 

5 

100 

DOT = DOT + A(I)*B(I) 

DOT 

6 


RETURN 

DOT 

7 

r 

END 

DOT 

8 


SUBROUTINE PLOAD(IDjFi>B»NNfP) 

PLOA 

1 


FORM LOAD UECTOR IN COMPACT FORM 

PLOA 

2 


DIMENSION ID(l)f F(1)»B(1) 

PLOA 

3 


DO 100 N=l> NN 

PLOA 

4 


J=ID(N) 

PLOA 

5 

100 

IF(J.GT.O) B ( J ) =F ( N ) -”P 

PLOA 

6 
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RETURN 

END 

C 

FUNCTION PRQPLD(T,J) 

PROPORTIONAL LOAD TABLE (ONE LOAD CARD ONLY) 

COMMON /CTDATH/ Q» HEAD(SO) , NUMNP, NUMEL, LAYER, NEQ, IPR 
DIMENSION ACS) 

IF <J ,LE. 0) GO TO 200 
C.,,. INPUT TABLE OF PROPORTIONAL LOADS 
1=1 

READC5, 1000) K,L,TMIN,TMAX, CA(KKK)»KKK=1»5) 

WRITECG, 2000) 0,HEAD, I,K,L,TMIN,TMAX, CA(KKK),KKK«1»5> 

RETURN 

C.... COMPUTE UALUE AT TINE T 
200 PROPLD =0.0 

IFCT.LT.TMIN ,OR, T.GT.TMAX) RETURN 
L « I1AX0CL, 1) 

PROPLD = Aa)+A(2)«T+A(3)w<SIN(A(4)»T'f-A(5)))*»L 
RETURN 

1000 FORNAT(HI5»?F10*0) 

2000 FORMAT ( A1 , 20A4//5X, 23HPR0P0RTI0NAL LOAD TABLE//11H NUMBER , 

1 43H TYPE EXP. MINIMUM TIME MAXIMUM TIME* 13X, 2HA1, 13X, 

2 2HA2, 1 3X, 2HA3 , 1 3X , 2HA4 , 1 3X , 2HA5/ ( 3 I 8 , ?G 1 5 , 5 ) ) 

END 

C 

SUBROUTINE PRTDISCUL, ID, X, D. F, T, NDM, NDF) 

C#»#» OUTPUT NODAL UALUES 
LOGICAL PCOMP 
COMMON /PROLuB/ i-’RQr 

COMMON /CTDATA/ Q» HEAD ( 30 )» NUMNP, NUMEL, LAYER, NEQ, IPR 
COMMON /LABELS/ PDIS(6),A(6),BC(2),DI(B),CDC3)»FD(3) 

COMMON /TMDATA/ TIME, DT, DDT, FORCE, ALPHA 

DIMENSION X(NDM, 1),B(1),UL(S) , IDCNDF, 1) ,F(NDF, 1)»T<1) 

DATA BL/4HBLAN/ 

DO 102 N=1 , NUMNP 
I F ( PCOMP ( X ( 1 , N ) , BL ) ) GO TO 101 
DO 100 1=1, NDF 
UL(I) = F(I,N)#PR0P 
K = IABSCID<I,M)> 

100 IF(K.GT.O) UL(I)=B(K) 

. T(N)=UL(3) 

101 CONTINUE 

102 CONTINUE 

WRITEC3,200l) (TCI), 1=1, NUMNP) 

RETURN 

2001 F0RMAT(GE12.4) 

END 

C 

SUBROUTINE FSTREA CUL, XL, LD, P, IX, ID, X, F, JDIAG, DR, D, NDF, NDM, NEN, 
a NST.NNEG) 

C**»* ELEMENT ROUTINE 

COMMON /CTDATA/ 0, HEAD (20), NUMNP, NUMEL, LAYER, NEQ, IPR 

COMMON /ELDATA/ N,NEL,MCT 

COMMON /ISHIDX/ ISW 

COMMON /PPDLOD/ PROP 

DIMENSION UL(NDF, 1), XL (NDM, 1) ,LD(NDF , 1),P( 1), IX (NEN, 1 ), 

1 ID (NDF, I ) , X(NDM, i),F(NDF, 1), JDIAG (i),DR(l)»B(l)»S(l) 

IF ( ISW.EQ.5) CALL PLOAD( ID, F, DR, NNEQ, PROP) 

MCT=Q 

DO 110 N=l, NUMEL 

CALL PFORM(UL,XL, LD, IX, ID, X, F, B, NDF, NDM, NEN, ISW) 

CALL ELMT01 (UL, XL, IXC 1 , N) , P, NDF, NDM, NST, ISW) 

IF(I5W.NE.4) CALL ADDSTF(DR,S, P, JDIAG, LD, 1,NEL*NDF, .FALSE, ) 

110 CONTINUE 
RETURN 
END 
C 

SUBROUTINE PFORMCUL, XL,LD, IX, ID, X, F. U, NDF, NDM, NEN, ISW) 

C*### FORM LOCAL ARRAYS 

COMMON /ELDATA/ N,NEL,MCT 
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PLOA 7 
PLOA 0 

PROP 1 
PROP r! 
PROP 3 
PROP 4 
PROP i) 
PROP tJ 
PROP ? 
PROP U 
PROP S 
PROP 10 
PROP 1* 
PROP lil 
PROP 13 
PROP 14 
PROP IS 
PROP IB 
PROP 17 
PROP 10 
PROP iii 
PROP CO 
PROP 21 

PRTD 1 
PRTD 2 
PRTD 3 
PRTD 4 
PRTD 5 
PRTD G 
PRTD 7 
PRTD 8 
PRTD 3 
PRTD 10 
PRTD 11 
PRTD 12 
PRTD 13 
PRTD 14 
PRTD 15 
PRTD IG 
PRTD IF 
PRTD 18 
PRTD 13 
PRTD 20 
PRTD 21 
PRTD 22 

FSTR 1 
FSTR 2 
FSTR 3 
FSTR 4 
FSTR 5 
FSTR G 
FSTR 7 
FSTR 8 
FSTR 3 
FSTR 10 
FSTR 11 
FSTR 12 
FSTR 13 
FSTR 14 
FSTR 15 
FSTR 1G 
FSTR IF 
FSTR ID 

PFUR 1 
PFOR 8 
PFOR 3 
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COMMON /PROLOD/ PROP 

DIMENSION ULCNDF, i ) » XLCNDM, 1 ) » LDCNDF, 1 ) * IXCNEN, 1 ) , IDCNDF, 1 ) , 
a XCNDM, 1)»F(NDF» 1)»U(1) 

DO 108 I-l.NEN 
II « IX(I»N) 

IF(II .NE. 0) GO TO 105 
DO 103 JaUNDM 

103 xlcj.i) » 0. 

DO 104 Jal,NDF 

ulcj*d a o. 

104 LUC Ji I) « 0 
GO TO ioa 

105 IID a II#NDF - NDF 
NEL * I 

DO 108 J*li NDM 
10G XLC J, I) « XCJ,II) 

DO 107 J=l» NDF 
K » lABSCIDCJ,!!)) 

ULCJ,I) a FCJ*II)#PROP 
IFCK.GT « 0) ULCJ, I)*U(K) 

IFCISN.EQ.G) K=IID+J 

107 LD(J,I) = K 

108 CONTINUE 
RETURN 
END 

C 

SUBROUTINE ELMT01 CUL, XL, IX, P, NDF, NDM, NST, ISM) 

C**»# LINEAR ELASTIC IN-PLANE a BENDING ELEMENT ROUTINE 
LOGICAL TAN 

COMMON /ELDATfl/ N,NEL,MCT 

COMMON /MTDATA/ RHO, UU12, El , E2? G12, G13* G23, THK, MIDTH 
COMMON /CQMP5T/ ABDCG,G),DS(2,2),QBR(3»3,25),QBS(2,2,25), 
a THC2S), 2KC25) 

COMMON /DMATIX/ DC10),DBCG,G),LINT 
COMMON /TNDATA/ TIME, DT, DDT, FORCE, ALPHA 
COMMON /GAUSSP/ SG<1G),TG(1G),MGC1S) 

COMMON /EXTRAS/ TAN 

DIMENSION ULCNDF, 1 ) , XL (NDM, 1 ) , IXC1 ) , PC 1 ) , SHPC3, 12), 

1 SIGTC3) , SIGBC3) ,SIGS(2), EPTC3) , EPBC3) , EPSC2) 

C 

DO 20 L^l, NST 
20 PCD = 0.0 

C.,.« COMPUTE NEUTRAL STRAINS AND STRESS RESULTANTS 
L a DCl) 

IFCISM.EQ.4) L=DC3) 

CALL PGAUSS(L,LINT) 

DO GOO Lai, LINT 

C .. COMPUTE ELEMENT SHAPE FUNCTIONS 

CALL SHAPECSGCL) , TGCL) , XL, SHP, XSJ, NDM, NEL, IX, .FALSE. ) 

C .. COMPUTE STRAINS AND COORDINATES 
DO 410 I«l,3 
EPTCI) a 0.0 
410 EPBCI) « 0.0 
DO 420 I«l, 2 
420 FPSCI) a 0,0 
XX a o.O 
YV a o.O 
DO 430 J«1»NEL 
XX a XX + SHPC3, J)»XL(1, J) 

YY a YY + SHPC3, J)«XLC2,J) 

C .. IN-PLANE STRAINS 

EPT(l) « EPT(l) + SHPC1,J)*UL(1.J) 

EPTC2) = EPTC2) + SHPC2, J)«UL(2, J) 

EPTC3) a EPTC3) + SHPC1, J)«UL(2» J) + SHPC2, J)#UL( 1, J) 

C .. BENDING CURUATURES 

EPBCI) a EPBCI) - SHPC1,J)#ULC4,J) 

EPBC2) a EPBC2) - SMPC2, J)»ULC5, J) 

EPBC3) a EPBC3) - SHPC1, J)*ULC5, J) - SHP(2, J)#UL(4, J) 

C .. SHEARING STRAINS 

EPS Cl) a EPS(l) + SHPC1, J)*UL(3» J) - SHPC3, J)*UL(4» J) 


PFOR 4 
PFOR 5 
PFOR G 
PFOR 7 
PFOR 8 
proR 8 
PFOR 10 
PFOR 11 
PFOR 12 
PFOR 13 
PFOR 14 
PFOR 15 
PFOR 18 
PFOR 17 
PFOR 18 
PFOR 19 
PFOR 20 
PFQR 21 
PFOR 22 
PFOR 23 
PFOR 24 
PFOR 25 
PFOR 2G 
PFOR 27 
PFOR 28 

ELMT 1 
ELMT 2 
ELMT 3 
ELMT 4 
ELMT 5 
ELMT G 
ELMT 7 
ELMT 8 
ELMT 9 
ELMT 10 
ELMT 11 
ELMT 12 
ELMT 13 
ELMT 14 
ELMT 15 
ELMT 1G 
ELMT 17 
ELMT 18 
ELMT 19 
ELMT 20 
ELMT 21 
ELMT 22 
ELMT 23 
ELNT 24 
ELMT 25 
ELMT 2G 
ELNT 27 
ELNT 20 
ELMT 29 
ELMT 30 
ELMT 31 
ELMT 32 
ELNT 33 
ELMT 34 
ELMT 35 
ELMT 3G 
ELMT 37 
ELMT 38 
ELMT 39 
ELMT 40 
ELMT 41 
ELMT 42 
ELMT 43 
ELMT 44 
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SHPC3»J)#ULC5»J) 

i»s) 


430 EPSC2) * EPSC2) + SHPC2,J)*ULC3,J) 
IFCISU*EQ*5.AND.TAN) 

* WRITEC9.90Q1) N,L. CEPBCII) »n«l,3>» CEPSCII), I 
3001 F0RMAT(2IB»5E12.4) 

: .. COMPUTE 3TRESS RESULTANTS 

DO 440 I a l * 3 
SIGT(I) “ 0, 

SIGBCI) » 0, 

DO 440 J»l»3 
SIGTCI) ■ SIGTCX) + 

440 SIGBCI) " SIGB(I) + 


ABDCI, J)#EPTCJ) + 
ABUCI+3»J)«EPTCJ) 


450 


DO 450 i~i,n 
SIGSCI) * 0. 

DO 450 J«l,2 

SIGSCI) a SIGSCI) + DSCI, J)#EPSCJ) 
IFCISW.GT.4) GO TO G20 

OUTPUT STRESS RESULTANTS AND STRAINS 


ABDC I, J+3)#EPBC J) 

+ ABDCI+3,J+3)WEPBCJ) 


MCT “ MCT - 2 


470 


470 


I l • * 
620 


IFCMCT.GT.O) GO TO 
WRITECG»2001) TIME 
MCT o 30 

WRITECB, 2002) N» XX. YY. EPT. EPB, EPS, SIGT» SIGB, SIGS 
GO TO BOO 

COMPUTE INTERAL FORCES 
DU n XSJ«WGCL) 


J1 

DO 


« 1 
G10 


PCJ1 ) 
PCJ1+1) 
pcjua) 

PCJl+3) 


J=1»NEL 


PCJ1 ) 
PCJl+l) 


PCJl+3) 

VcJl+4) * PCJ1+4) 


CSHPC 1, J)*SIGT(1 )+SHPC2, J)#5I,GTC3) )*BU 
CSHPC2, J)#SIGTC2)+SHPC 1, J) W SIGTC3) )»BU 
CSHPC 1. J)*SIG3C 1 )+3MPC2, J)*3IG3C2) )*DU 
CSHPC 1, J)#SIGB( 1 )+SHP(2» J) W SIGB(3)+SHP(3» J) 
*SIGSC1))»DU 

CSHPC2,J)«SIGBC2)+SHPC1,J)**SIGBC3)+SHPC3,J) 

«SIGSC2))#DU 


BIO J1 » J1 + NDF 
SOO CONTINUE 
RETURN 
C 

2001 FQRMATC1H1// 

^ SX.SHTIME °» E12 »3//5X, 33HELEMENT STRAINS/STRESS RESULTANTS// 

1 8H ELEMENT , 3X, 7H1-CQQRD, 3X» 7H2-C00RD, 4X, SHXX-STRAIN. 4X» 

2 9HYY-STRAIN>4X»9HXY"STRAIN» 3X» 1QHKXX-STRAIN.3X, 

3 10HKYY-STRAIN, 3X» 10HKXY-STRAIN, 4X, SHSX-STRAIN, 4X, 

4 9HSY-STRAIN/28X.BCGX.7H-STRESS)/) 

2002 FORMAT C I0» 2F10.4, 8E13.4/28X.8E13.4) 

END 

C 

SUBROUTINE PGAUSSCLL.LINT) 

C«»** GAUSSIAN POINTS AND WEIGHTS FOR TWO DIMENSIONS 
COMMON /GAUSSP/ SGClG).TGaG).WG(lG) 

DIMENSION LR(9),LZCS),LWC9),WR(2),GRC2),GCC2) 

DATA LR/-i,l,l,-l»0,l»0»-l,0/,L2/-l,-l,l»i,-l,Q,l«0,Q/ 

DATA LW/4<*25, 4*40,64/ 

DATA GR/O. 8S113G31 1594053, 0.339981043584S5G/ 

DATA GC/1. 0,0, 3333333333/ 

DATA WR/O , 347854845137454, 0 . G521451548G254G/ 

LIMT = LL*LL 
L«IABSCLL) 

GO TO Cl, 2, 3, 4), L 
C.... 1X1 INTEGRATION 

1 SGC 1 ) « 0. 

TGC 1 ) « 0. 

WGC1) ■ 4. 

RETURN 

C. . . , 2X2 INTEGRATION 


2 G » 1./SQRTC3.) 
IFCLL.LT.O) G~l» 
DO 21 1=1,4 
SGC I) » G*LRCI) 
TGC I) ■ G«L2(I) 


ELMT 45 
ELMT 40 
ELMT 47 
ELMT 415 
ELMT 49 
F.tNT 
ELMT bl 
ELMT liC 
EUH U3 
ELMT 1* 
ELMT UJ 
CLf ir \r 
ELMT tii’ 
ELMT H) 
ELMT ;»a 
ELMT Mil 
ELMT OS 
ELMT UP 

elmt r:n 

ELMT Gi 
ELMT DM 

am L'i 

ELMT 07 
elmt nn 

ELMT (551 
ELMT 70 
ELMT ?1 
ELMT ?H 
ELM f ?3 
ELMT 74 
ELMT 75 
ELMT ?G 
ELMT 7r‘ 
ELMT 7B 
ELMT 79 
ELMT 80 
ELMT 01 
ELMT 82 
ELMT 03 
ELMT 84 
ELMT 03 
ELMT QG 
ELMT 07 
ELMT 83 
ELMT 83 
ELMT 90 

PGAU 1 
PGAU 2 
PGAU 3 
PGAU 4 
PGAU 5 
PGAU G 
PGAU 7 
PGAU 8 
PGAU 9 
PGAU 10 
PGAU 11 
PGAU 12 
PGAU 13 
PGAU 14 
PGAU 15 
PGAU 1G 
PGAU 17 
PGAU 10 
PGAU 19 
PGAU 20 
PCAU 21 
PGAU 22 
PGAU 23 
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21 UGCI) * 1. 

RETURN 

C.... 3X3 INTEGRATION 

3 G = SORTCO.G) 

IFCLL.LT.O) G-l. 

H = l./Ol. 

DO 31 1*1,9 
SG( I ) * G*LRCI) 

TGCI) » G*LZCI) 

31 HG(I) » 

RETURN 

C.... -1X4 INTEGRATION 

4 DO 41 1=1, 4 

11 n 1+M0DCI+1.2) 

12 * 1 

IFCI.GT.2) 12 - 2 
DO 41 J=l»4 
JJ « (I-l)*4+J 
SGCJJ) * LR(J)«GRCI1) 

IFCLL.LT.O) SGCJJ) « LRCJ)*GCCI1) 

TGC JJ) = LZC J)#GR( 12) 

IFCLL.LT, 0) TGCJJ) » LZ(J)#GCCI2) 

41 WGCJJ) * WR(I1)«WR<I2) 

RETURN 

END 

C 

SUBROUTINE SHAPECS5, TT, X, SHP, XSJ, NDM, NEL, IX,FLG) 

C«*#» SHAPE FUNCTION ROUTINE FOR TWD DIMENSIONAL ELEMENTS 
LOGICAL FLG 

DIMENSION SHPC3,4)»XCNdM* i)»§C4)»TC4)»X3C2*2),8XC2»2),IHCS) 

DATA S/-0.5»0.5,Q.5,-0.5/,T-'-0.5,-0.5,0.5,0.5/ 

C.,.. FORM 4-NODE QUADRILATERIAL SHAPE FUNCTIONS 
DO 100 1=1,4 

SHPC3, 1 ) = ( 0 . S+S ( I ) **SS ) * ( 0 . 5+T ( I ) *TT ) 

SHPC1, I) = SCI)«C0.5+TCI)*TT) 

100 SHPC2, I ) = T(I)*C0.5+SCI)#SS) 

IFCNEL.GE.4) GO TO 120 

C.... FORM TRIANGLE BV ADDING THIRD AND FOURTH TOGETHER 
DO 110 1=1,3 

110 SHP(I»3) = SHP(I,3)+SHP(I,4) 

C.... ADD QUADRATIC TERMS IF NECESSARY 

120 IFCNEL.GT.4 .AND. NEL.LT. 10) CALL SHAP2CSS,TT,SHP, IX,NEL) 

C.... ADD CUBIC TERMS IF NECESSARY 

IFCNEL.GT.S) CALL SHAP3CSS, TT, SHP, IX, NEL) 

C.... CONSTRUCT JACOBIAN AND ITS INUERSE 
DO 130 1=1, NDM 
DO 130 J=l, 2 
XSCI,J) = 0.0 
DO 130 K=1»NEL 

130 XS(I,J) = XS(I,J)+ XCJ*K)*SHPCI»K) 

XSJ = XSC 1, 1 )^XS(2,2)-XSC1, H)«XS(2, 1) 

IFCXSJ ,GT. 0.00000001) GO TO 135 
URITEC6, 2000) IX 
C STOP 

135 IFCFLG) RETURN 

SXC1,1) = XS(2,2)/XSJ 
SX<2,2) = XSCl.lVXSJ 
SXC1,2) = -XSC1,2)/XSJ 
SXC2, 1) = -XSC2, 1)/XSJ 
C.... FORM GLOBAL DERIUATIUES 
DO 140 1=1, NEL 

TP = SHP(l,I)wSXCl,l)+SHP(2,I)«SXC2,l) 

SHP (2, 1 ) = SHPa,I)#SXCl,2)+SHP(2,msX<2,2) 

140 SHPC 1, I ) = TP 
RETURN 

2000 FORMAT C5X, G7H««SHAPE ERROR 01««- ZERO OR NEGATIUE JACOBIAN DET. FOR 
-'ELEMENT NODES: /20X, i2I4) 

END 

C 

SUBROUTINE SHAP2CS,T,SHP, IX, NEL) 
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origins pri ; 

OF POOR 


C»»** ADD QUADRATIC FUNCTIONS AS NECESSARY 

SHAP B 

DIMENSION IX(9),SHP(3» IE) 

SHAP 3 

SS - (l.“S#3)/2» 

SHAP 4 

T2 * Cl.-T*T)/2. 

SHAP 5 

DO 100 I«5> NEL 

SHAP 0 

DO 100 J*l* 3 

SHAP 7 

100 SHPCJ, I) » 0.0 

SHAP 0 

C,... MIDSIDE NODES (SERENDIPITY) 

SHAP 9 

IF(IX(5).EQ.O) GO TO 101 

SHAP 10 

SHP(1,S) « -S*C1.-T> 

SHAP 11 

SHP(2,5) * -SB 

SHAP 12 

SHPC3.5) - S2*K1.-T) 

SHAP 13 

101 IF(NEL.LT.G) GO TO 107 

SHAP 14 

IFdX(G).EQ.O) GO TO 102 

SHAP 15 

SHP(l.G) - T2 

SHAP 1G 

SHP(2,6) » ~T«(1,+S) 

SHAP 17 

SHP(3» G) - T2»(l.+S) 

SHAP 18 

102 IFCNEL.LT. 7) GO TO 107 

SHAP 19 

IF(IX(7) .EQ.Q) GO TO 103 

SHAP EO 

SHP(1»7) - -S*(l.+T> 

SHAP 21 

SHP(B»7) « se 

SHAP VP, 

SIIP (3. 7) = 82*(1.+T) 

SHAP 23 

103 IF(NEL.LT.G) GO TO 107 

SHAP 84 

IF(IX(8).EQ.O) GO TO 104 

SHAP 85 

SHPC1,8) « -T2 

SHAP 2S 

SHP(2»8) « -T«(l.-S) 

SHAP 27 

SHP(3, 8) - T2»C1.-S) 

SHAP 28 

C.... INTERIOR NODE (LAGRANGIAN) 

SHAP 29 

104 IF(NEL.LT.S) GO TO 107 

SHAP 30 

IF(IXCS).EQ.O) GO TO 107 

SHAP 31 

SHP Cl, 9) = -4.*S»T2 

SHAP 32 

SHPC2.3) « -4.*T*52 

SHAP 33 

SHP(3,3) a 4»*S2#T2 

SHAP 34 

C. . . . CORRECT EDGE NODES FOR INTERIOR NODE (LAGRANGIAN) 

SHAP 35 

DO 10G J“l» 3 

SHAP 3G 

DO 105 I®1.4 

SHAP 37 

105 SHP(JiI) = SHPCJ, I) - 0.25«SHPCJ,9) 

SHAP 38 

DO 10G 1*5.8 

SHAP 39 

10G IF(IXd).NE.O) SHP(J.I) = SHP(J.I) -0.5*SHPC J. 3) 

SHAP 40 

C.... CORRECT CORNER NODES FOR PRESENCE OF MIDSIDE NODES 

SHAP 41 

107 K = 8 

SHAP 42 

DO 103 1=1.4 

SHAP 43 

L = I + 4 

SHAP 44 

DO 108 J=l» 3 

SHAP 45 

108 SHP(J.I) = SHP(J.I) - 0.5#(SHP(J.K)+SHP(J»L) ) 

SHAP 4G 

109 K = L 

SHAP 47 

RETURN 

SHAP 48 

END 

SHAP 4S 

c 

SUBROUTINE SHAP3(S.T, SHP. IX. NEL) 

SHAP 1 

C«**« ADD CUBIC FUNCTION AS NECESSARY (SERENDIPITY) 

SHAP 2 

DIMENSION IX(12).SHP(3. 12) 

SHAP 3 

DO 100 1=5, NEL 

SHAP 4 

DO 100 J=l» 3 

SHAP 5 

100 SHP< J, I)=0. 0 

SHAP G 

IFdX(5) .EQ.O) GO TO 101 

SHAP 7 

Sl=-l./3. 

SHAP 8 

Tl=-1. 

SHAP 9 

CALL CSHAPE(S,T, Sl.Tl.SHP, 1,5) 

SHAP 10 

101 IFdX(G).EQ.O) GO TO 102 

SHAP 11 

Sl=l. 

SHAP 12 

Tl=-l./3. 

SHAP 13 

CALL CSHAPECS. T» Si , T1 , SHP, 2. G) 

SHAP 14 

102 IF(IX(7) .EQ.O) GO TO 103 

SHAP 15 

Sl*l./3» 

SHAP J.B 

Tl=l. 

SHAP 17 

CALL CSHAPE(S»T, Sl.Tl.SHP, 1,7) 

SHAP 18 

103 IFdX(8).EQ.O) GO TO 104 

SHAP 13 

Sl=-1. 

SHAP 20 

Tl=l./3. 

SHAP 21 
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CALL C5HAPECS, T * 51 , T1 » SHP, 2, 8) 

104 IFCIXC9) .EQ.O) GO TO 105 
S1--1, 

T1&-1./3. 

CALL CGHAPECS, T» 51 , T1 » SHP, 2*9) 

105 IFCNEL.LT . 10) GO TO 500 
IFCIX(IO).EQ.O) GO TO 10G 
Sl=l./3. 

Tl**-1. 

CALL CSHAPECS»T»S1»T1»SHP» 1,10) 

106 IPCNEL.LT .11) GO TO 200 
IFCIXCID.EQ.O) GO TO 107 
Sl*l. 

Tl=l,/3, 

CALL CSHAPECS»T*S1»T1»5HP»E» 11) 

107 IFCNEL.LT. 12) GO TO 200 
IFCIXC 12) .EQ. 0). GO TO 200 
S1--1./3. 

Tl«l. 

CALL CSHAPECS,T»S1,T1»SHP» 1, 12) 

C.... CORRECT CORNER NODES 
200 DO 210 1=1.4 
11=1+4 
12=1+8 

IFCI.EQ.l) 13=1+7 
IFCI.GT.l) 13=1+3 
IFCI.LT.4) 14=1+9 
IFCI.EQ.4) 14=1+5 
DO 210 d=i, 3 

210 SHPCJ,I)=SHPCJ,I)-2./3.*CSHPCJ,ID+SHPCJ,i£))-l./3,#CSHPCJ,I3> 
~ +SHPCJ, 14) ) 

RETURN 

END 

C 

SUBROUTINE CSHAPECS, T, 51, Tl, SHP, K, L) 

C*«*« SUPPLEMENTAL ROUTINE FOR THE SHAPE FUNCTIONS 
DIMENSION SHPC3.12) 

C=9./32. 

GO TO C1,2),K 

1 SHP C 1 » L ) =C# ( 1 . +T1*T ) * ( 9. *Sl-2 . *5-27 . *S1*S*S ) 
SHPC2,L)=C*T1#C1.-S«S)«C1.+9.#S1*S) 

SHP C 3 , L ) =C# C 1 . +T1*T ) * Cl . -S«S ) # C 1 . +9 . *S1 *S ) 

RETURN 

2 SHPCl»L)=C*Sl*Cl.-TwT)«a.+9,*Ti*T> 
SHPC2,L>=DKl,+Sl#S)#C9.«Tl-2.tfT~27.*Tl#T*T) 
SHPC3,L)=C«Cl.+SlwS)»(l,-T^T)wCl.+9.*TH*T) 

RETURN 

END 

C 

SUBROUTINE PMESHCIDL, XL, IX, ID, X, F, JDIAG, NDF, NDM, NEN, NKM) 

C*«*n INPUT MESH DATA 
LOGICAL PRT»ERR»PCOMP 

COMMON /CTDATA/ 0, HEADC20 ) , NUMNP, NUMEL, LAYER, NEQ, IPR 
COMMON /MTDATA/ RHO, UU12, El, E2,G12,G13,G23,THK, WIDTH 
COMMON /LABELS/ PDISCG),ACG),BCC2),DICS),CDC3),FDC3) 

COMMON /EXDATA/ 0LAWC4) 

COMMON /RODATA/ UR,IQ,NDS 

DIMENSION IDL(G),XLC7), IXCNEN, 1), IDCNDF, 1),XCNDM, 1), 
a F ( NDF , 1 ) , DUM C 1 ) > WDC 13) » JDI AG C 1 ) 

DATA WD/4HCOOR, 4HELEM, 4HMATE, 4HBOUN, 4HF0RC* 4HROD , 

^ 4HEND , 4HPRIN, 4HNOPR, 4HPAGE, 4HEXPE/ 

DATA BL/4HBLAN/, LIST/11/, PRT/. TRUE. / 

C.... INITIALIZE ARRAYS 
ERR = .FALSE. 

DO 501 1=1,4 
501 OLAWC I )=0. 

DO 502 N=l, NUMNP 
DO 502 1=1, NDF 
IDCI,N)=0 
FCI,N)=0. 
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CSHA 9 
CSHA 10 
CSHA 11 
CSHA 12 
CSHA 13 
CSHA 14 

PMES 1 
PMES 2 
PMES 3 
PMES 4 
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nr 

OF POOR QUALITY 

HO 

Ii 

502 CONTINUE 

PMES 22 

rn.it 

tu 

! j; 

C. . . . READ A CARD AND COMPARE WITH MACRO LIST 

PMES 23 

»» 

10 READC5, 1000) CC 

PMES 24 


DO 20 1-1, LIST 

PMES 85 


20 IF(PCOMP<CC* MD( I ) ) ) GO TO 30 

Piles 08 

ill 

GO TO 10 

PMES 8? 

i ?! 

30 GO TQ (1,2, 3, 4,5, 6,7,0,9,11,12), I 

PMES <B 


C.... NODAL COORDINATE DATA INPUT 

PMES OS 


1 DO 102 N“1»NUMNP 

PMES 30 


102 XC1,N)« BL 

Fiu'S 31 

• '• 

CALL GENUEC(NDM» XL, X, CD, PRT » ERR) 

PMES 88 

isv 

GO TO 10 

PMES 33 


C,..* ELEMENT DATA INPUT 

PMES 31 

mm 

2 L*0 

PMES 33 


DO 20G I=1*NUMEL»50 

PMES 36 

v*& 

IF(PRT) WRXTE(6>2001) 0,HEAD, <K*K=1>NEN) 

PMES 37 

J * MINOCNUMEL.I+49) 

PHrT) ,5U 


DO 20G N=I» J 

FHSS 3S 

mm 

: '> 

IF(L-N) 200,202,203 

PMES 40 

,| 

200 READC5, 1001) L» (IDL(K),K=1,MEM)»LX 

PMES 41 


IFCL.F-Q.O) L-NUMEL+1 

PMES 42 


IF(LX.EQ.O) LX*1 

FM3 <’3 

«**•* 

IF(L-N) 201,202,203 

PMES 44 

■ H 

201 WRITE(6»3001 ) L,N 

PMES 43 

i i 

mu 

ERR * .TRUE. 

PMES 4G 

GO TQ 20G 

PMES 4? 


202 NX « LX 

PMES 43 

mtn 

DO 20? K“1»NEN 

PMES 49 

h* 

207 IX(K»L) * IDL(K) 

PMES 50 

v L 

GO TO 205 

503 IX<NEn»n) = IX(NEN.N-l) 

PMES 51 


PMES 52 

rr*. 

DO 204 K=1»NEN 

PMES S3 

7 f: 

IX(K,M) * IX(K,M~1) + NX 

PMES 04 

' i 1' 

St a 

204 IF(IX(K»N-l).EQ.O) IX(K,H? * 0 

PMES 53 

205 IF(PRT) WRITECG, 2002) N, ClX(K,N),K=i,NEM) 

PMES 5G 


20G CONTINUE 

PMES 5? 

j 

■|l 

GO TO 10 

PMES S3 


C. . . . MATERIAL DATA INPUT 

PMES 59 


3 WRITE(6»2004) O.HEAD 

PMES GO 


CALL MATLIB 

PMES 61 

. i. 

GO TO 10 

PMES G2 

Ilf 

C.... READ IN THE RESTRAINT CONDITIONS FOR EACH NODE 

FMES G3 

» ii 

4 IF(PRT) URITE(6, 2000) 0* HEAD, (I,BC» 1=1, NDF) 

PMES 64 


N * 0 

PMES G5 


NG = 0 

PMES GG 

<!p 

420 L = N 

PMES G? 

ill' 

i& it 

LG = NG 

PMES G8 

READ(5» 1001 > N,NG,IDL 

PMES G9 


IF(N.LE. 0 .OR. N.GT.NUMNPj GO TO 50 

PMES 70 

mm 

DO 41 1=1, NDF 

PMES 71 

■ If 

ID(I»N) = IDL(I) 

PMES 72 

T» kt 

41 IFCL.NE. 0 .AND, IDLCD.EQ.O .AND. IDCI,L) .LT.Q) IDCI,N)=-1 

PMES 73 


LG = ISIGN(LGiN-L) 

PMES 74 


42 L - L+LG 

PMES 75 

■ j r 

IF((N“L)«LG .LE. 0) GO TO 420 

PMES 7G 

* J. 
(11: 

teM 

DO 43 1=1, NDF 

PMES 77 

43 IF(ID(I»L-LG) .LT. 0) ID(I,L) = -1 

PMES 78 


GO TO 42 

PMES 79 

mm 

50 DO 48 N=1,NUMNP 

PMES 80 

A 

DO 4G 1=1, NDF 

PMES 81 


46 IF(ID(I*N) .NE, 0) GO TO 47 

PMES 82 


GO TO 48 

PMES 83 

! '{ ft 

4? IF (PRT) WRITECG, 2007) N, (ID(I»N)» I=1»NDF) 

PMES 84 

48 CONTINUE 

PMES 85 

h 

GO TO 10 

PMES 03 

C.... FORCE/DISPL DATA INPUT 

PMES 87 


5 CALL GENUEC ( NDF , XL , F , FD , PRT , ERR ) 

PMES 03 

\ > ii 

GO TO 10 

PMES 89 

ill 

C. . . . END OF MESH DATA INPUT 

PMES 90 

J-s 

C COMPUTE THE PROFILE OF GLOBLE ARRAYS 

PMES 91 

n 

oi 

m 


[if 


7 IP* (ERR) STOP 
CAUL PROFILUBIAG,ID.XX*NDF,NEN*NKM*PRT> 

RETURN 

Cm,. PRINT OPTION 
0 PRT n .true. 

GO TO 10 

Cm., NOPRINT OPTION 
3 PRT » .FALSE. 

GO TO 10 

Cmm READ IN PAPER EJECTION OPTION 

11 READ(5#i000) 0 
GO TO 10 

Cm., INPUT EXPERIMENTAL INDENTATION LAW 

12 READC5. 1007) (QLflU(X)t I«l*4) 

HRITE<6.20Q8) OrHEAD. (QLAW(I)* I»l*4) 

GO TO 10 

C.,m INPUT INITIAL IMPACT CONDITION 
(5 HR1TE(G»2009) 0*HEAD 
REAQ(5» 1002) NQ* INDF*UR 
WRITER* 8010) NO* INDF.UR 
F(INDF.NQ)«1,Q 
IQalD(INDF.NQ) 

GO TO 10 

C.... INPUT/OUTPUT FORMATS 

1000 FORMAT < A4 * 7SX * A 1 ) 

1001 FORMAT? iGXB) 

1002 FORMAT (SIS. F10.0) 

1007 FORMAT? 4F1Q.0) 

2000 FORMAT?Al»20A4//SX»lOHNOfiAL B.C«*7X//8X»3HN0D& »9(X7>A4. A2)/1X) 

2001 FORMAT cnif 20 A4//SX * 8HELEMENT9//3X» r HELEMENT* 

a 14(13. SH M0I)E)/(2QX* 14Q3.5H NODE))) 

aooa FORMATdiOi 1,4I8/'(10X* 1418)) 

8004 FQRMAT?A1»20A4//SX* IMMATERIAL PROPERTIES) 

2007 FORMAT? I10» 9X13) 

2003 FORMAT?A1.20A4//6X» ^EXPERIMENTAL INDENTATION LAW*// 

1 1QX* /CONTACT COEFFICIENT: **E12,4/ 

2 1 OX* /CRITICAL INDENTATION: P EXE. 4/ 

3 1 OX .^CONSTANT S: N E12.4/ 

4 10X.NPOWER INDEX OF UNLOADING LAN: N F12.3) 

3001 FORMAT (5X * B6H«*»PME3H ERROR 01*w ELEMENT* IS* 

a SSI-1 APPEARS AFTER ELEMENT* 15) 

2003 F0RMAT(A1*20A4*//5X* /IMPACT OF LAMINATED PLATEN) 

2010 PORMAT(//10X* /IMPACT NODAL POINT* *110/ 

a 10X.NIMPACT D.O.F.* NiIlO/ 

a 10X, /INITIAL IMPACT UEL0CITV:N*Eia.4) 

END 

C 

SUBROUTINE GENUCCCNDM* XL* X. CD. PRT* ERR) 

GENERATE: REAL DATA ARRAYS BY LINEAR INTERPOLATION 
LOGICAL PRT » ERR* PCONP 

COMMON /CTDATfV 0* HEAD ( 20 ), NUMNP* NUMEU LAYER* NEQ* XPR 
DIMENSION X(NDM»1)*XL(7)*CD(3) 

DATA BL/4HBLAN/ 

N»0 

HG'-O 

102 L«N 
LG*NO 

READ(S.IOOO) N.NG*XL 

IF(N*LE.0 .OR, N.GT. NUMNP) GO TO 108 

DO 103 1*1, NON 

103 X(I*N)«XL(I) 

IF(LG) 104,10a»104 

104 LG“ISIGN(LG»N~L) 

LI»(IABS(N-L+LG)-1)/IABS(LG) 

DO 105 1*1, NDH 

105 XL(I)*?XCI,N)-X?I,L>>/U 
10G L“L+LG 

IF((N-L)«LG .LE. 0) GO TO 103 
IFCL.LE. 0 .OR. L.GT. NUMNP) GO TO 110 
DO 107 I n l. NDM 
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PMES 32 
PMES 93 
PMES 34 
PMES 95 
PMES 98 
PMES 97 
PMES 99 
PMES 99 
PME5100 
PMES101 
PMESlOa 
PMES103 
PMES104 
PMES 105 
PMES10S 
PMES107 
PMES1QB 
PMES109 
PMES 110 
PI1ES111 
PMES112 
PMES113 
PMES114 
PMES 115 
PMES11G 
PMES 117 
PMES118 
PMES119 
PMES120 
PMES1S1 
PMESiaa 
PMES123 
PMES 124 
PMES1S5 
PMES12B 
PMES137 
PMESia8 
PMES129 
PMES 130 
PMES 131 
PMES132 
PMES, 133 
PMES134 
PMES 135 
PMES13G 
PMES137 

GENU 1 
GENU a 
GENU 3 
GENU 4 
GENU 5 
GENU 6 
GENU 7 
GENU 8 
GENU 9 
GENU 10 
GENU 11 
GENU IS 
GENU 13 
GENU 14 
GENU 15 
GENU 16 
GENU 17 
GENU 18 
GENU 19 
GENU ao 
GENU 21 
GENU 22 
GENU 23 
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107 Xa.L)*Xa,L-LG)+XLCI) 

GO TO 10G 

110 WRITE(B,3Q00) L, (CDCI), 1*1,3) 

ERR « .TRUE. 

GO TO 10E 

10B DO 109 I*1*NUNNP»50 

IFCPRT) WRITE CG, 2000)0, HEAD, <CDCL),L«1,3), CL,CDC1)»CDC2),L*X.NDM> 

M * MIND CNUMNP* 1+49) 

DO 1C3 J*I,N 

IF<PCOMPCXa,J),BL) .AND. PRT) URITECG,2008) N 
109 IP C . NOT . PCOtIP C X ( 1 » J ) » BL ) . AND . PRT > WRITE<6,2009) J* CX(U J),L*1»NBM) 
RETURN 

1000 FQRMAT(2I5* 7P10 , 0) 

8000 FORMAT < A1 » 20A4,V5X* 5HN0DflL.3fl4//BX»4HN0DE»9(X7»fl4t A2)) 

2Cj8 FORMAT C5X,21H»«GENUEC WARNING 01**, I 10, 

* 3SH HAS NOT SEEN INPUT OR GENERATED) 

2009 FQRMATC 110, 9F13.4) 

3000 FORMAT (5X»44H*«GENUEC ERROR 01»*ATTEMPT TO GENETATE NODE* 15, 

1 3H IN,3A4) 

END 

C 

SUBROUTINE PRQFILt JDIAG, ID* XXvNDF,NEN.NKM»PRT) 

C«»»* COMPUTE PROFILE OF GLOBAL ARRAYS 
LOGICAL PRT 

COMMON /CTDATA7 0, HEAD C 20 )*NUHNP,NUMEL, LAYER, NEQ* XPR 
DIMENSION JDIAGC 1 ) , IDCNDF, 1 ) , IX(NEN, 1 ) , EQC2) 

DATA EQ/4H DOF,2H. ✓ 

C.... SET UP THE EQUATION NUMBERS 
NEQ « 0 

DO 50 N*1,NUNNP 
DO 40 1*1, NDF 
J * IDCItN) 

IFCJ) 30,20,30 
20 NEQ » NEQ + 1 
IDCItN) ■ NEQ 
JDIAGCNEQ) * 0 
GO TO 40 
30 ID(I,N) * 0 
40 CONTINUE 
50 CONTINUE 

IF C. NOT. PRT) GO TO 70 
WRITE(G,2000> 0»HEAD, (I,EQ, 1*1, NDF) 

DO GO I«1,NUMW» 

GO WRITEC6* 2001 ) I, CIDCK, t),K“l,NDF) 

C.... COMPUTE COLUMN HEIGHTS 
70 DO 500 N=1»NUNEL 
DO 400 1*1, NEN 
II « IX(I,N) 

IF C 1 1 .EQ, 0) GO TO 400 
DO 300 K*1»NBF 
KK * IDCK, II) 

IFCKK.EQ. 0) GO TO 300 
DO 200 J=I,NEN 
JJ * IXtJ.N) 

IFCJJ.EQ.O) GO TO 200 
DO 100 L=l, NDF 
LL * ID(L, JJ) 

IFCLL.EQ.O) GQ TO 100 
M * MAX0(KK,LL) 

JDIAGCM) * NAXOC JDIAG(M), IAB5CKK-LL)) 

100 CONTINUE 
200 CONTINUE 
300 CONTINUE 
400 CONTINUE 
500 CONTINUE 

C.... COMPUTE DIAGONAL POINTERS FOR PROFILE 
NKM « 1 
JDIAGC 1) * 1 
IFCNEQ.EQ.l) RETURN 
DO GOO N=E» NEQ 
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GENU 24 
GENU Kfi 
GENU Eli 
GENU Dr' 
GENU 28 
GENU pn 
GENU 30 
GENU 31 
GENU I'D 
GENU 33 
GENU 34 
GENU 33 
GENU 3U 
GENU 37 
GENU 30 
GENU 30 
GENU 40 
GENU 41 
GENU 42 
GENU 40 

PRO! 1 
PROF n 
PROF 0 
PROF 4 
PROF S 
PROF G 
PROF ? 
PROF 0 
PROF 8 
PROF 10 
PROF 11 
PROF ID 
PROF 13 
PROF 14 
PROF 15 
PROF IS 
PROF 1? 
PROF 18 
PROF 13 
PROF 20 
PROF 21 
PROF as 
PROF 23 
PROF 24 
PROF Do 
PROF SG 
PROF 27 
PROF DU 
PROF 29 
PROF 30 
PROF 31 
PROF 32 
PROF 33 
PROF 34 
PROF 35 
PROF 3G 
PROF 37 
PROF 38 
PROF 39 
PROF 40 
PROF 41 
PROF 42 
PROF 43 
PROF 44 
PROF 45 
PROF 46 
PROF 47 
PROF 48 
PROF 49 
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COMP.GHOSITE / 


GOO JDIAGCN) * JDIAGCN) + JDIAGCN-l) + 1 
NKM ■ JDIAGCNEQ) 

5000 FORNATtftl*BOM//3X*lBHEQUftTXON NUMBERS//GX* 5HN0DE * 
a 9<XS»A4*fl8)/lX) 
aooi FORMAT(iio»aui) 

RETURN 
EN1) 

C 

SUBROUTINE MTLXB 

MATERIAL PROPERTIES ROUTINE 
COMMON /CTUATfV 0» HEADC80) * NUMNP* NUMCL* LAYER* NEQ. I PR 
COMMON /NTDATA/ RHO»UUia» El *EB* GIB* G13*GB3»THK» WIDTH 
COMMON /COMPST/ ABDCGrG).DSCa*a).QBRC3,3,E5),QBSCa*e»a5)* 
a TH(85)*2KCa5) 

COMMON /DMATIX/ DU0).D8<G»B>*LXNT 
DIMENSION I-ID(S) 

DATA WD/6H ISO-.GH ORTHO* SHTR0PIC.6H 
C*... INPUT MATERIAL PROPERTIES 
READ(5* 1000) L1»LB>K.THK» WIDTH 
READCS* 1001 ) RH0»UU12»EliE5iG15iG13»GB3 
DO 150 J*l,3 
DO 150 1*1.3 

IFCI.EQ.3 .OR. J.EQ.3) GO TO 150 
DSC J* I) « 0. 

150 ABDCJ. I ) b ADDCJ+3.I) * ABDCJ. 1+3) * ABDCJ+3.I+3) - 0, 

LI « MIN0(4*MAX0C1»L1>> 

DC!) « LI 

L2 « NM0C4.NAX0a.L8)) 

DCS) a LB 
DC3) « K 
LINTaO 

iFCEi-ES) iao*no»i 80 
uo Gia«Ei/ca»»u.+uma)) 
ji-i $ ua«3 
go to aoo 
120 Jla4 $ ja^s 

XFCLAYER.EQ.l) Jl«2 $ J8“3 

S00 WR1TECG* 8000) LAYER. MDCJ1 ) * WDC JB) . THK. Ei» EH. GIB* G13* GB3»UU1S* 
a RH0*L1»LB» M . 

CALL CNPD 
RETURN 

C*m. FORMAT FOR INPUT-OUTPUT 

1000 FORMAT ( 3 IS* 8F1 0 * 0 ) 

1001 PORNATCrTlO.O) 

2000 FORMATC/5X* IB* 1BH LAYER (S) OF.aAG.BlH PLATE WITH THICKNESS* 

1 F10 . 4*V10X» ISHYOUMGKS MODULUS* 10X* XEiat*. E10 . 4. 10X. *EB»»X* E10 . 4/ 

2 10X* 15HSHEAR MODULUS* SX* ^G 1 B®M* E10 . 4* SX* ^G.13*!*** E10 . 4* SX* 

3 XG23"F*E10»4/10X* 15HP0ISS0N RATIO* BX.HUUlfiaX*FS.3/iOX* 

4 FHDENSITV* i?X.*RHOai»!,EiO,4/iOX. 13HGAUSS PTS/DIR. lEX.^Li**. IS. 

5 SX»FLBbH»IS/ 10X»1BHSTRESS POINT , 14X* s^K*^. 15/) 

END 

C 

SUBROUTINE CMPD 

C*«w* COMPUTE s^ABD* MATRIX AMD *DS;* MATRIX 

COMMON /CTDATA/ 0 * HEAD (80)* NUMNP* NUMEL* LAYER. NEQ. IPR 
COMMON /MTDATA/ RHQ.UU1B* El *EB*G18*G13#GB3» THK* WIDTH 
COMMON /COMPST/ ABDC6»6)»DS(8»8)*QBR(3*3»B5)»QBS<&»B»S5)i 
a THCaS). 2 K(aS) 

DIMENSION QC3*3),GBCa,B).TK(85) 

LLaLAYER 

NM-LL+l 

REWKUflOOO) (L.THCD.TKCL). I«1*LL) 

2R(l)«TTKaO.O 
DO IS IbI»lL 
TTIv^TTR+TK C I ) 

2K(I+l)aTKtt)*2Ka) 

15 CONTINUE 
DO *85 I u l* MM 
2K ( I ) »2K ( I ) -TTK/ 8 . 
as CONTINUE 
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PROF 50 
PROF 51 
PROF 58 
PROF 53 
PROF 54 
PROF 55 
PROF 5G 


NATL 2 
NATL 3 
NATL 4 
NATL 5 
MATL G 
NATL ? 
MATL 0 
MATL 9 
MATL 10 
MATL 11 
MATL 18 
NATL 13 
MATL 14 
MATL 15 
MATL IB 
MATL IF 
MATL 10 
MATL 19 
MATL SO 
NATL El 
NATL aa 
MATL 83 
MATL 84 
MATL 85 
MATL EG 
MATL a? 
MATL HO 
MATL 89 
MATL 30 
MATL 31 
MATL 33 
MATL 33 
MATL 34 
NATL 35 
MATL 36 
MATL 3? 
NATL 38 
MATL 39 
MATL 40 
MATL 41 
MATL 48 
MATL 43 

CMPD 1 
CMPD a 
CMPD 3 
CNPD 4 
CMPD 5 
CMPD 6 
CMPD 7 
CMPD 0 
CMPD 9 
CMPD 10 
CMPD a 
CMPD IB 
CMPD 13 
CMPD 14 
CMPD IS 
CMPD IS 
CMPD 1? 
CMPD 18 
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ORKSBW. P- 5 -®,};’ 
OF POOvl Qli&Ls n ^ 


DEL=4.*ATANC1.>/1S0. 

DEN ■ 1. - EG#UU12»*2/E1 

0(1*1) * El/DEN 

0(2,2) * E27DEN 

0(1,2) = 0(2,1) * VU12«Q<2,2) 

0(3,3) = G12 

Q(l» 3) « 0(2,3) » Q(3, 1) « 0(3,2) * 0.0 

0S(1,JL) = G13 

QS(2,2) * G23 

0S(1,2) * GS(2,1) * 0.0 

DO 40 1=1, LL 

ANGL-TH«)*DEL 

C=CQS(ANGL) 

W=SIN(ANGL) 

QBRd,l,I)=Qd»l)*C#M+2,#CQd,2)+2.*4K3,3)mC*W)**2-H3(2,2)*W**4 
QBRCl, 2, I)=QBR(2» l*I)=(Qd»l)+Q(2,2)-4»*Q(3»3) )tt(C«W) { t#2 
$ +Q<1»2)«<M#»4 +Cix*4 ) 

QBRC2, 2, I ) =0(1, l)tH|»«4+2.«(Q(l*2)+2.#Q(3*3) )*(C«W)fc»2+Q(2*2)«C#*4 
QBR(1»3, I)*QBR(3»l*I)«<G<l,l)-Q(l,2)-2.»Q(3,3))#M»OH»3 + 

$ (0(l,2)-0(2,2)-f'2.«0(3,3))4*C«Wi^3 

0BR(2,3,I)=QI3R(3,2,I)=(Q(l,i)*-Q(l,2)-2.^Q(3,3))«W»tt3 *»C+ 

$ (Q( 1, 2)~Q(2, 2H2,«0(3, 3) )#M«C«m»3 

QBR(3,3,I) = <Q(l.l)+QCa,2)-2.wQd,2)“2.-«Q(3,3))*KW»C)»tt2+ 

$ Q(3,3)*(W«”4 +C««4 ) 

QBS(1, 1,1) = QS(1,1)#C*h> 2 + QS(2,2)WRww2 
QBS(2,2, 1) = DS(1,1)«W»«2 + QS(2,2)»C*#2 
QB3(1»2, I) = QBS(B, 1,1) = (QS(1, 1>-QS<2,2))*C*W 
40 CONTINUE 
DO 50 J=l* 3 
DO 50 !<=1» 3 
DO 50 1=1, LL 

ABD( J ,K )= ABfl(J ,K )+QBR(J,K,I)*(2K(I+l)“ZKCI)) 

ABD( J+3, K )= ABD(J ,K+3)= ABD(J+3,!<)+QBR(J,K, I)* 

$ (ZK(I+l)»*2-ZK(I)*w2)/2. 

ABD( J+3, K+3)= ABD ( J+3 , K+3 ) +OBR ( J » K * I ) * ( ZK ( 1+1 ) w«3-ZK ( I ) w*3 ) 73 . 

50 CONTINUE 
DO 55 1=1, B 
DO 55 J=l» 6 

IFCI.GE.3 .OR, J.GE.3) GO TO 55 
IF(ABS(DS(I* J) ) .LT. 1.E-0G) DS(I,J)=0.0 

55 IF(ABS(ABD(I, J)) .LT. 1.E-0G) ABD(I,J)=0.0 
URITE(G,2001) ((ABD(I, J), J=l,6), 1=1, G) 

DO GO J=l, 2 
DO GO K=l, 2 
DO GO 1=1, LL 

GO DS( J,l<) = DS( J, K) + QBS(J»K, l)#(ZK(I+i)“ZK(l)) 

MRITE(6, 2002) ((DS(I, J), J=l,2), 1=1,2) 

1000 FORMAT(I5,F5.0,F10.0) 

2001 FORMAT (7/, IX, 10HABD MATRlX77G(2X, 6E13.47) ) 

2002 F0RMAT(/,1X,9HDS MATRIX772(2X, 2E13.47) ) 

RETURN 

END 

C 

SUBROUTINE KMLIB 
C*hh*« ASSEMBLE GLOBLE ARRAY 
COMMON Gd) 

DIMENSION Md) 

EQUIUALENCE (G(1),M(1)) 

COMMON 7IS14IDX7 ISU 

COMMON /CTDATA7 0, HEAD(20 ) , NUMNP, NUMEL, LAYER, NEQ, IPR 
COMMON /LODATA/ NDF * NDM, NEN, NST , NKM 
COMMON /PARATS7 NPAR(14) » NEND 
N1=NEND 

N2=N1+NST«NST*IPR 
IF( IGU.LE.2) NE=N2+NKM*IPR 
IF(ISW.GT»2) NE=N2+NEQ*IPR 
CALL SETMEM(NE) 

CALL PZERO ( G ( NEND ) , NE-NEND ) 

CALL MASS01(G(NPAR(1) ),G(NPAR(2)),M(NPAR(3)),G(NPAR(4) ), 

I M(NPAR(5)),M(NPAR(G)),G(NPAR(7))»G(NPAR(8)),M(NPAR(9)), 


CMPD 19 
CMPD 20 
CMPD 21 
CMPD 22 
CMPD 23 
CMPD 24 
CMPD 25 
CMPD 25 
CMPD 2? 
CMPD 28 
CMPD W 
CMPD 30 
CMPD *J1 
CMPD 32 
CMPD 33 
CMPD 34 
CMPD 35 
CMPD 3S 
CMPD 37 
CMPD 39 
CMPD ;;s 
CMPD 40 
CMPD 41 
CMPD 42 
CMPD 43 
CMPD 44 
CMPD 45 
CMPD 4G 
CMPD 47 
CMPD 40 
CMPD 43 
CMPD 50 
CMPD 51 
CMPD 52 
CMPD 53 
CMPD 54 
CMPD 55 
CMPD 5S 
CMPD 57 
CMPD 58 
CMPD 59 
CMPD GO 
CMPD SI 
CMPD G2 
CMPD S3 
CMPD G4 
CMPD G5 
CMPD GG 
CMPD G7 
CMPD G8 
CMPD G9 
CMPD 70 

KMLI 1 
KMLI 2 
KMLI 3 
KMLI 4 
KMLI 5 
KMLI G 
KMLI 7 
KMLI 8 
KMLI 9 
KMLI 10 
KMLI 11 
KMLI 12 
KMLI 13 
KMLI 14 
KMLI 15 
KMLI IB 
KMLI 17 
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2 GCNPAR(ll) )>GCN1)*GCN2),NDF,NDM»NEN,N5T»NKM) 

RETURN 

END 

C 

SUBROUTINE MAS5Q1CUL, XL,LD» P, IX, ID, X, F, JDIAG, B, S, A, NDF, NDN, NEN, 
a NST,NKM) 

C**#» FORM MASS MATRIX 

COMMON /CTDATA^ 0, HEAD (20 ) , NUMNP, NUMEL, LAVER, NEQ, IPR 

COMMON /MTDATA7 RHO, UU12, El , E2» G12, G13, G23, THK, WIDTH 

COMMON /DMATIX7 DC10),DBC6,G),LINT 

COMMON /ELDATA/ N,NEL,MCT 

COMMON /ISWIDX/ ISM 

COMMON /GAUSSP/ SG(16),TGC16),WGC16) 

DIMENSION ULC 1 ) , XLCNDM, 1 ) , LDCNDF, 1 ) , PC 1 ) . IXCNEN, 1 ) , IDCNDF, 1 ) , 

1 XCNDM, I ) ,FC1 ) » JDIAGC 1 ) ,BCi) »SCNST , 1), ACl) »SHP(3, 12) 

C.... LOOP ON ELEMENTS 
DO HO N«l, NUMEL 
DO 10 1=1, NST 
DO 10 J=1,NST 
10 SCI,J)=0. 

C.... SET UP LOCAL ARRAYS 

CALL PFORMCUL, XL, LD, IX, ID, X, F, B, NDF, NDM, NEN, ISM) 

C.... COMPUTE CONSISTENT MASS MATRIX 
L = DC 1 ) 

CALL PGAU5SCL, LINT) 

DO 500 L=1,LINT 

C .. COMPUTE SHAPE FUNCTIONS 

CALL SHAPE C SG C L ) , TG ( L ) , XL , SHP , XS J , NDM , NEL , IX * . FALSE . ) 

DU b WG(L)*XSJ*RHO*THK 

C .. FOR EACH NODE J COMPUTE DB=RHO»SHAPE*DU 
K1 = 1 

BO 500 J=1,NEL 
Mil = SHPC3, J)*DU 
W33 = Hll*THK»*2/12. 

C .. FOR EACH NODE K COMPUTE MASS MATRIX CUPPER TRIANGULAR PART) 
J1 = K1 

DO 510 K=J,NEL 

SCJ1 ,K1 ) = SCJ1 ,K1 ) + SHPC3, K)#W11 

SC Jl+3,Kl+3) = SC Jl+3,Kl+3) + SHPC3,K)»W33 
510 J1 = J1 + NDF 
500 K1 = I<1 + NDF 

C .. COMPUTE MISSING PARTS AMD LOWER PART BY SYMMETRY 
NSL = NEL*NDF 
DO 530 K=1,NSL»NDF 
DO 520 J=K,NSL,NDF 
SC J+2, K+2) = SC J+1,K+1) = SCJ »K ) 

SC J+4, K+4) = SC J+3,K+3) 

SCK ,J ) = SCJ , K ) 

SCK+3, J+3) = SCJ+3,K+3) 

SC K+2, J+2) = SCK+l, J+l) = SCJ ,K ) 

520 SCK+4, J+4) = SC J+3, K+3) 

530 CONTINUE 

IFCISW.EQ.2) GO TO 100 
C.... LUMPED MASS MATRIX 
SUM 1 =0.0 
SUM2 =0.0 
SUMD1 =0.0 
SUMD2 =0.0 
DO 540 1=1, NSL, NDF 
SUMD1 = SUMD1 + SCI, I) 

SUMD2 = SUMD2 + SC 1+3, 1+3) 

DO 540 J=1,NSL,NDF 

sum = sum + sa,j) 

540 SUM2 = SUMS + SC 1+3, J+3) 

DO 550 1=1, NSL, NDF 
PCI) = S(I,I)«SUM1/SUMD1 
PCI+2) = PCI+1) = PCI) 

PCI+3) = SCI+3,I+3)*SUM2/SUMD2 
550 PC 1+4) = PCI+3) 

C.... ADD TO TOTAL ARRAY 


KMLI 18 
KMLI 19 
KMLI 20 

MASS 1 
MASS 2 
MASS 3 
MASS 4 
MASS 5 
MASS G 
MASS 7 
MASS 8 
MASS 9 
MASS 10 
MASS 11 
MASS 12 
MASS 13 
MASS 14 
MASS 15 
MASS 1G 
MASS 17 
MASS 18 
MASS IS 
MASS 20 
MASS 21 
MASS 22 
MASS 23 
MASS 24 
MASS 25 
MASS 2G 
MASS 27 
'MASS 28 
MASS 29 
MASS 30 
MASS 31 
MASS 32 
MASS 33 
MASS 34 
MASS 35 
MASS 36 
MASS 37 
MASS 38 
MASS 39 
MASS 40 
MASS 41 
MASS 42 
MASS 43 
MASS 44 
MASS 45 
MASS 46 
MASS 47 
MASS 48 
MASS 49 
MASS 50 
MASS 51 
MASS 52 
MASS 53 
MASS 54 
MASS 55 
MASS 56 
MASS 57 
MASS 58 
MASS 59 
MASS GO 
MASS 61 
MASS 62 
MASS G3 
MASS 64 
MASS 65 
MASS 66 
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100 CALL ADDSTF(A.S.P.JDIAG»LD.NSTrNEL*NDF». FALSE.) 

110 CONTINUE 
REWIND 2 

IFCISW.EQ.2) WRITEC2) CACI), I*1»NKM) 

IF( ISW.EQ.3) WRITEC2) CACI), 1*1, NEQ) 

RETURN 

END 

C 

SUBROUTINE RDDIPCT 

C«*»# 

LOGICAL FLAG 
COMMON G(l) 

DIMENSION MCI) 

EQUIUALENCE (GCD.MCl)) 

COMNON /CTDATA/ 0» HEAD C 20 ) » NUMNP, NUMEL, LAYER# NEQ, I PR 

COMMON /LODATA/ NDF,NDM»NEN»NST» NKM 

COMMON /PARATS/ NPAR(14)»NEND 

COMMON /RODATA/ UR.IQ.NDS 

COMMON /ROELEM/ NER,NEQR,ER 

DATA FLAG/ . FALSE . A NER/20/, ER/30000000 . / 

IF (FLAG) GO TO 50 
NEQR=2#(NER+1) 

NKMR=7»NER+3 

N1=NEND 

N2=N1+NEQ*IPR 

N3=N2+NEQ»IPR 

N4=N3+NEQ«IPR 

N5=N4+NKMR^IPR 

NG'-NS+NEORwIPR 

N7=NG+NEQR 

NB=N7+NEQR»IPR 

M9=N8+NEQR#IPR 

N10=N9+NEQR«IPR 

N11=N10+NEQR#IPR 

NE=M11+NEQR»IPR 

CALL SETMEMCNE) 

CALL PZERO ( G ( NEND ) » NE-NEND ) 

FLAG*. TRUE. 

50 CALL WIHPCT(G(NPAR(l))»G(NPAR(2))»M(NPAR(3))»G(NPARC4))f 

1 MCNPARC5)), MCNPARCG) )» GCNPARC7)), GCNPARC8) )» 

2 H ( NP AR ( 9 ) ) » G ( NP AR C 1 0 ) ) » G ( NP AR ( 1 1 ) ) » G ( N 1 ) , G ( N2) » 

3 GCN3),G(N4),GCN5),M(N6),GCN7),G(N8),GCN9),GCN10), 

4 GCN11)) 

RETURN 

END 

C 

SUBROUTINE WIMPCT(UL»XL,LD, P, IX, ID,X,F» JDIAG,DR,U» B» U* A.RK.RM. 
a JDR.RU.RU.RA.RB.FR) 

C»*** SOLUE IMPACT PROBLEM 
LOGICAL FLAG. TAN 
COMMON GC1) 

DIMENSION MCI) 

EQUIUALENCE CGCD.MCl)) 

COMMON /CTDATA/ 0, HEADC20 ) . NUMNP. NUMEL. LAYER* NEQ. IPR 

COMMON /TMDATA/ TIME, DT. DDT. FORCE, ALPHA 

COMMON /LODATA/ NDF, NDM, NEN, N5T » NKM 

COMMON /NITERS/ ITR 

COMMON /PARATS/ NPARC14),NEND 

COMMON /RODATA/ UR.IO.NDS 

COMMON /ROELEM/ NER.NEQR.ER 

COMMON /CONSTS/ A0.A2.A4. A5, AG, A7» A8, AREA 

COMMON /PROLOD/ PROP 

COMMON /ISWIDX/ ISW 

COMMON /EXTRAS/ TAN 

DIMENSION ULC 1 ) » XLC 1 ) , LDC 1 ) , P C 1 ) , IXC 1 ) , IDC1 ) » XC 1 ) » FC 1 ) » JDIAG Cl). 

1 DRCD.UCD.BCl), UC1), ACD.RKCDtRMCl), JDRCD.RUCl), 

2 RUC1)»RAC1),RBC1)»FR(1)»QC3),QPC3) 

DATA ITR/5/, FLAG/. FALSE. /, HIL/1. 4/, INTE/24/ 

IF (FLAG) GO TO 50 

DO 1 1=1,3 


MASS G7 
MASS 63 
MASS G3 
MASS 70 
MASS 71 
MASS 72 
MASS 73 

RODI 1 
RODI 2 
RODI 3 
RODI 4 
RODI S 
RODI 6 
RODI 7 
RODI 8 
RODI 0 
RODI 10 
RODI U 
RODS 12 
RODI 13 
RODI 14 
RODI 15 
RODI IS 
RODI 17 
RODI 13 
RODI 19 
RODI 20 
RODI 21 
RODI 22 
RODI 23 
RODI 24 
RODI 25 
RODI 26 
RODI 27 
RODI 28 
RODI 29 
RODI 30 
RODI 31 
RODI 32 
RODI 33 
RODI 34 
RODI 35 
RODI 3G 
RODI 37 

WIMP 1 
WIMP 2 
WIMP 3 
WIMP 4 
WIMP 5 
WIMP 6 
WIMP 7 
WIMP 8 
WIMP 9 
WIMP 10 
WIMP 11 
WIMP 12 
WIMP 13 
WIMP 14 
WIMP 15 
WIMP 1G 
WIMP 17 
WIMP 18 
WIMP 19 
WIMP 20 
WIMP 21 
WIMP 22 
WIMP 23 
WIMP 24 
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0(13=0.0 

WIMP 25 


QPCI)°0.0 

WIMP 2B 

1 

CONTINUE 

WIMP 2? 


IU5-1 

WIMP 28 


TAN=. FALSE. 

WIMP 29 


REWIND 2 

WIMP 30 


READC23 (B(I),I=1»NEQ) 

WIMP 31 


FGRCEaO.O 

WIMP 32 


ALPHAS. 0 

WIMP 33 


PROP=0. 0 

WIMP 34 


NNEQ=NDF*NUMNP 

WIMP 35 


A0=S./'<WIL*DT)##2 

WIMP 3G 


A2=6./CWIL*DT) 

WIMP 37 


A4=A0/WIL 

WIMP 38 


A5=-A2/WIL 

WIMP 39 


AB=1 .-3./WIL 

WIMP 40 


A7=DT/2. 

WIMP 41 


A8=DDT/G. 

WIMP 42 


CALL FORMROD ( RK » RM , JDR ) 

WIMP 43 


DO 10 I=l,NEOR 

WIMP 44 

10 

RUCI3=-UR 

WIMP 45 


QC2)=-UR 

WIMP 46 


FLAG=.TRUE. 

WIMP 47 

50 

ISW=5 

WIMP 48 


IFCIDS.EO.NDS) TAN=.TRUE. 

WIMP 49 


CALL FSTREA ( UL, XL, LD» P, IX, ID, X, F» JDIAG, DR, U, NDF, NDM* NEN, NST, NNEO) 

WIMP 50 


DO 20 1=1, NEQ 

WIMP 51 


A«)=DR(I)/B<I> 

WIMP 52 


UCI)=UCI)+DT*ACI) 

WIMP 53 


UCI)=U(I)+DT*UCI) 

WIMP 54 

20 

CONTINUE 

WIMP 55 


QPC1)=UCIQ) 

WIMP 5B 


0PC2)=U(IQ) 

WIMP 57 


GPC3)=ACI0) 

WIMP 58 


DO 30 I=1,NEQR 

WIMP 59 


RB ( I 3 =RM (I ) « ( A0#RU ( I 3 + A2*RU ( I ) +2 . *RA (I 3 3 

WIMP BO 

30 

CONTINUE 

WIMP 61 


RB I Q=RU ( 1) +DT*RU C 1 3 +DDT/3 . «RA C 1 3 

WIMP 62 


ROT=0. 000001 

WIMP 63 


ICOU=0 

WIMP G4 


DO 100 IT=1» ITR 

WIMP G5 


RUT=RBIQ+Q(3)»DDT/ / B. 

WIMP BG 


AF=-RUT-QPC13 

WIMP G7 


CALL RODLOADCFIQ, AF3 

WIMP G8 


DO 110 I=1,NEQR 

WIMP 69 


FRCI3=RBCI3 

WIMP 70 

110 

CONTINUE 

WIMP 71 


FR ( 1 3 =FR < 1 3 + C 1 . - W I L 3 *FORCE+ W I L*FI Q 

WIMP 72 


CALL ACTCOLCRK.FR, JDR, NEQR, .FALSE. , .TRUE. , 0) 

WIMP 73 


Q C 3 ) =A4« C FR ( .U -RU ( 1 ) ) +A5«RU C 1 ) +A6-»RA C 1 3 

WIMP 74 


RUTT=RBIQ+QC33*DDT/6. 

WIMP 75 


ROTR=ABS ( ( RUTT-RUT 3 /RUTT 3 

WIMP 7B 


IFCROTR.LT.ROT) ICOU=l 

WIMP 77 


IFCIC0U.GT.03 GO TO 200 

WIMP 78 

100 

CONTINUE 

WIMP 79 

200 

DO 210 1=1, NEQR 

WIMP 80 


FR ( I ) =A4« ( FR C 1 3 -RU CD) +A5*RU ( I ) +AG*RA C 1 3 

WIMP 81 


RUCI3=RU(I3 +DTWRU C I ) +AS« C FR ( I ) +2 . «^RA CD) 

WIMP 82 


RU ( I ) =RU C 1 3 +A7« C FR C 1 3 +RA CD) 

WIMP 83 


RACI)=FRCI) 

WIMP 84 

210 

CONTINUE 

WIMP 85 


oci)=rucd 

WIMP SB 


Q C 2 ) =RU Cl) 

WIMP 87 


0(33=RAC 1 ) 

WIMP 88 


FORCE=FIQ 

WIMP 89 


PROP=FORCE 

WIMP 90 


ALPHA=-0 Cl) -OP Cl) 

WIMP 91 


RODFR=RU CINTE3 *AREA*ER 

WIMP 92 


WRITEC8, 8001 ) FORCE, ALPHA, RODFR, COC I ) , 1=1, 3) 

WIMP 93 

8001 

FORMAT C BE 12. 4) 

WIMP 94 
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WIMP 35 


IF(IDS.GT.NDS) IDS-1 

WIMP 96 


TAN». FALSE. 

WIMP 97 


RETURN 

WIMP 30 

p 

END 

HIMP 93 

w 

SUBROUTINE F0RMR0D(RK» RM» JDR) 

form y 

o*** 

FORM STIFFNESS AND MASS MATRICES OF ROD 

ror;si s 


COMMON /RODATA/ UR.IQ.NDS 

FORM J 


COMMON /ROELEM/ NER.NEQR.ER 

FORM 4 


COMMON /CONSTS/ AQ» A2» A4, A5» AS. A?» A8» AREA 

FORM U 


DIMENSION RKC1),RM(1),JDR(2),D(G) 

FORM t; 


DATA RHQR7. 00032257, RL/1.07 

FORM 7 


DATA D/.22, .30. .43. .48, .50. .G25/ 

FORM B 


EL-RL/HER 

FORM 9 


PAI=4. #ATAN(1. ) 

FORM 10 


JDR< 0=1 

FORM 11 


JDR(2)=3 

FORM 12 


DO 100 1=1. NER 

FORM 13 


IF(I.LT.B) A=PAI#(D(I)/2.)#*2 

FORM 14 


IF(I.GE.B) A=PAI«(D(G)/2.)«*2 

FORM IS 


TT=A*ER/30./EL 

FORM lb 


Jl=2»(I+l)~l 

FORM l? 


J2=J1+1 

FORM 10 


J1M1=J1~1 

FORM 19 


J1N2«J1-B 

FORM 20 


JDR(Jl)=JDR(JlMl)+3 

FORM SI 


JDRC J2)=JDR< J1 )+4 

FORM PS 


Ki s JBR( J1M2) 

FORM 03 


K2=JDRCJ1N1)~1 

FORM SI 


RK(I<1 )=RK(K1 )+TT*3B. 

FORM S3 


RK(K2 )=RK(K2 )+TT*3.«EL 

FORM CO 


RK ( K2+1 )«RK ( KH+1 ) +TTM . *EL**2 

FORM 27 


RK ( K2+2 ) “RK ( K2+H ) -TT *3B . 

FORM SB 


RK ( K2+3 ) = RK < K2+3 ) -TT»3 . <*EL 

FORM 83 


RK ( K2+4 ) “RK < K2+4 ) +TT”36 . 

FORM 30 


RK ( K2+5 ) =RK ( K2+5 ) +TT#3 . *EL 

FORM 31 


RK ( K2+B ) =RK < K2+B ) -1 T*EL*#2 

FORM 32 


RK ( K2+7 ) =RI< ( K2+7 ) -TT*3 . «EL 

FORM 33 


RK ( K2+0 )=RK< KH+8 ) +TT*4 . «EL»*2 

FORM 34 


TT=RHQR*A*EL 

FORM 33 


L1=2*I-1 

FORM 3G 


RMCL1 )=RM(L1 )+TT/2. 

FORM 37 


RM(Ll+l>=RM(Ll+n+TT»EL*#2/420. 

FORM 30 


RM(Ll+2)=RM(Ll+2)+TT/2. 

FORM 39 


RM(Ll+3)=RM(Ll+3)+TT«EL»*2/420. 

FORM 40 

100 

CONTINUE 

FORM 41 


AREA=A 

FORM 42 


DO 20 1=1. NEQR 

FORM 43 


J=JDR(I) 

FORM 44 

20 

RK(J)=RK(J)+AO#RM(I) 

FORM 45 


CALL ACTCOL(RK. RM. JDR, NEQR, . TRUE. , .FALSE. ,0) 

FORM 46 


RETURN 

FORM 47 

r 

END 

FORM 48 

w 

SUBROUTINE RODLOAD (F.AF) 

RODL 1 

C**#* 

COMPUTE CONTACT LOADING 

RODL S 


LOGICAL RELD.UNLD.PIL 

RODL 3 


COMMON 7TMDATA7 TIME, DT, DDT, FORCE, ALPHA 

RODL 4 


COMMON /EXDATA/ 0(4) 

RODL 5 


DATA UNLD/ . FALSE . PIL/ . FALSE . /, RELD/ . FALSE. / 

RODL 6 


IFCPIL) GO TO 10 

RODL 7 


AMAX=AMIN=FMAX=0 . 0 

RODL 8 


PIL=.TRUE. 

RODL 9 

10 

IF (RELD) GO TO 50 

RODL 10 


IF(UNLD) GO TO 20 

RODL 11 


F=Q(1)*AF*«1.5 

RODL 12 


IF(F.GE. FORCE) RETURN 

RODL 13 


UNLD=.TRUE. 

RODL 14 


AMAX=ALPHA 

RODL 15 
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FMAX=FORCE 

IF0WlX.GT.QC2)) UK*FMAX/ C C 1 . -Q C 3 ) ) # AMAX+Q < 2 ) *Q C 3 ) )**Q C 4 ) 
I F C AMAX . LE . Q C 2 ) ) UK»FMAX/AMAX*#Q C 4 ) 

AMIN=Q(3)«-(AMAX-QC2) ) 

IFCAMIM.LT. 0. ) AMIN=*0.0 
20 IFCAF.LE.AMIN) GO TO 30 
F*UK<KAF-AMIN)««QC4) 

IFCF.LT. FORCE) RETURN 
RELD-.TRUE. 

RK=FHAX/ C AH AX-AN IN ) 1 . 5 
50 IFCAF.LE.AMIN) GO TO 30 
F=RK« ( AF- AN I N ) «wi . 5 
RETURN 
30 F=0.0 
RETURN 
END 


RODL 1G 
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RODL 2G 
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RODL 29 
RODL 30 
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